Finding accreted stars in the Milky Way: clues from NIHAO simulations

Exploring the marks left by galactic accretion in the Milky Way helps us understand how our Galaxy was formed. However, finding and studying accreted stars and the galaxies they came from has been challenging. This study uses a simulation from the NIHAO project, which now includes a wider range of chemical compositions, to find better ways to spot these accreted stars. By comparing our findings with data from the GALAH spectroscopic survey, we confirm that the observationally established diagnostics of [Al/Fe] vs. [Mg/Mn] also show a separation of in-situ and accreted stars in the simulation, but stars from different accretion events tend to overlap in this plane even without observational uncertainties. Looking at the relationship between stellar age and linear or logarithmic abundances, such as [Fe/H], we can clearly separate different groups of these stars if the uncertainties in their chemical makeup are less than 0.15 dex and less than 20% for their ages. This method shows promise for studying the history of the Milky Way and other galaxies. Our work highlights how important it is to have accurate measurements of stellar ages and chemical content. It also shows how simulations can help us understand the complex process of galaxies merging and suggest how these events might relate to the differences we see between our Galaxy's thin and thick disk stars. This study provides a way to compare theoretical models with real observations, opening new paths for research in both our own Galaxy and beyond.


INTRODUCTION
The history of the Milky Way is a puzzle that has taunted astronomers for decades.The long lifetimes and rather conserved chemical composition of stars makes them key pieces in this puzzle, allowing us to use their star light as fossil record to gain insights into the historical processes that have led to our Galaxy as we know it today (Freeman & Bland-Hawthorn 2002).The advent of new large-scale stellar surveys, providing a wealth of information including measurements of elemental abundances 1 for millions of stars, has enabled insights into the long and complex history of our Galaxy that have been previously unattainable (Jofré et al. 2019).Astrometric information for more than 1.5 billion stars from the Gaia satellite (Gaia Collaboration et al. 2016Collaboration et al. , 2018Collaboration et al. , 2021) ) has given us a new understanding of the structure of our Galaxy, while spectroscopic surveys like the Galactic Archaeology with HERMES (GALAH) Survey (da Silva et al. 2015) or the Apache Point Observatory Galactic Evolution Experiment (APOGEE, Majewski et al. 2016) are complementing this vision with chemical fingerprints for millions of stars.Most notably, these surveys have led to the confirmation of a major accretion event logarithmic ratios of two elemental number densities  X and  Y compared to the Sun, that is, [X/Y] = log 10 . around 8-10 billions years ago by Belokurov et al. (2018) and Helmi et al. (2018).While this major merger of the Gaia-Sausage-Enceladus (GSE) and the early Milky Way could explain the bimodality of stars in the Milky Way's disk, the causality is yet to be established beyond doubt.Our best avenue is to quantify the chemical and dynamical properties of accreted stars to answer outstanding questions like how much gas the GSE merger brought in to explain the significantly different chemistry of young versus old disk stars.
From our vantage point within the disk-dominated region of the Milky Way, it is difficult to identify the relatively low number of accreted stars.Intensive research has tested a variety of ways to find such stars -an extensive list of which is presented by Buder et al. (2022).A first impulse was given by the series of papers by Nissen & Schuster (2010, 2011); Nissen et al. (2014) and Schuster et al. (2012) who analysed the chemistry of stars in the kinematic halo of the Milky Way.They identified a population within the Galactic halo that showed significantly lower enhancements across a number of elements, hinting to an extra-Galactic origin.Starting from dynamic arguments, Belokurov et al. (2018) and Helmi et al. (2018) rediscovered this population as stars with outstanding eccentric and radial orbits and proved its significance via dynamical simulations and chemical evolution arguments.Especially the inference of merger masses via chemical patterns is only possible because of our ability to observe surviving dwarf galaxies (Tolstoy et al. 2009;Helmi et al. 2018;Hasselquist et al. 2021;Carrillo et al. 2022).Peaking further into the halo, Naidu et al. (2020) then showed that the halo is entirely comprised of substructure, including accreted stars.Studies of extragalactic sources can now also spatially resolve these source broadly and use age-metallicity relations to identify accreted stars (Martig et al. 2021), with tens of edge-on galaxies being observed by the MUSE integral field spectrograph through large programs like GECKOS (van de Sande et al. 2023).Such analyses are however limited by the fact that spatial and dynamical selections are less robust in a spatially and dynamically evolving galaxy over billions of years.
The key advantage of chemical tagging, when compared to dynamical tagging via velocities or energies, is that the composition of heavy elements, once locked in the atmosphere of a star, is expected to be much more conserved in a dynamically evolving galaxy.Building upon this premise, Hawkins et al. (2015), Das et al. (2020) and Buder et al. (2022) assessed a number of potential abundance combinations for chemically tagging accreted stars.Among the many elements in the periodic table, they identified especially Na, Al, and Cu as exceptionally useful, if they can be measured.Using these elements together with the purer tracers of core-collapse supernovae (CCSN), like Mg, and supernovae type Ia (SNIa), like Mn, they found the abundance planes of [Na/Fe] vs. [Mg/Mn] and [Al/Fe] vs. [Mg/Mn] to reveal accreted stars in rather isolated loci.
Ages and abundances are, however, not directly observable and difficult to obtain from stellar spectra and any age-inference tool.Significant progress has been made in terms of the number and quality of observations, but we are still limited by our ability to observe a significant number of accreted stars with accurate and precise ages and abundances.While the GALAH survey can for example measure Na abundances, it is particularly difficult to measure Al abundances from its spectra.The observed stellar ages of GALAH DR3 are even more uncertain -with an average age uncertainty of 33 ± 20%.Cosmological simulations, such as those used in this study, carry no age uncertainties and the particular simulation (Buck et al. 2021) we use here, introduced in Sec.2.2, now also includes a wide range of chemical abundance information.
Cosmological simulation projects such as FIRE (Bonaca et al. 2017;Horta et al. 2023), EAGLE (Mackereth et al. 2019), NIHAO-UHD (Buck 2020;Buck et al. 2021Buck et al. , 2023)), AURIGA (Grand et al. 2020), VINTERGATAN (Agertz et al. 2021) or HESTIA (Khoperskov et al. 2023b) allow us to follow a self-consistent growth and chemical enrichment history, which is particularly important for metal-poor stars formed at redshifts  > 1. Simulation groups have made further huge progress in reaching the numbers of stellar particles necessary to split datasets in phase-space and abundance space.The latest versions of cosmological simulations thus provide an exceptional data set for exploring the changing chemodynamic environment of our Milky Way -which is neither in equilibrium (see e.g.Antoja et al. 2018;Bland-Hawthorn et al. 2019) nor had an isolated evolution (Helmi 2020).We can thus study the theoretically expected impact of massive mergers on the chemical evolution of Milky Way-like galaxies (Buck et al. 2021(Buck et al. , 2023) ) as well as our observational diagnostic tools without uncertainty margins to consider.
The aim of this project is to assess diagnostic tools for identifying accreted structures from observations in the Milky Way by comparing them with theoretical predictions from cosmological simulations.This assessment will be based on our ability to observe and quantify clear patterns and trends.In addition to testing chemical tagging, that is, relying only on abundances, we also test the use of age-abundance planes and make predictions on the necessary uncertainties for the diagnostic properties when assuming the Milky Way analogue as ground truth.
The paper is structured as follows: Sec. 2 describes the data of GALAH observations and NIHAO simulations.Sec. 3 performs an initial comparison of diagnostic plots, using different abundance planes.Sec. 4 considers how stellar abundances of accreted and insitu stars change over time, and quantifies the results.Sec. 5 discusses our findings and their implications.Sec.6 concludes the paper.

Observational data from the GALAH Survey
Observations are provided through the GALactic Archaeology with HERMES (GALAH) survey.For our study, we use the abundance measurements from the third data release (DR3) of the GALAH Survey with its 588,571 stars (Buder et al. 2021).The main program of the survey targets nearby stars via their visual magnitudes (typically 9 <  < 12 and 12 <  < 14) while neglecting the Galactic plane (|| > 10 deg).However, the data release also includes auxiliary programs of the K2 and TESS footprints (Sharma et al. 2018(Sharma et al. , 2019) ) as well as open and globular clusters (for more details see Buder et al. 2021).Stars were targeted with the 2dF fibre optics positioner system at the Anglo-Australian Telescope (Heĳmans et al. 2012;Farrell et al. 2014) and the light of up to 400 stars was simultaneously fed into the High Efficiency and Resolution Multi-Element Spectrograph (HERMES, Barden et al. 2010;Sheinis et al. 2015).HERMES delivers high-resolution ( ∼ 28, 000) spectra in four wavelength bands across the optical range which were reduced with the pipeline by (Kos et al. 2017).Stellar parameters ( eff , log , [Fe/H],  mic ,  broad , and  rad ) and abundances for up to 312 different elements across multiple nucleosynthesis channels are estimated using a modified version of the spectrum synthesis code Spectroscopy Made Easy (sme Valenti & Piskunov 1996;Piskunov & Valenti 2017) and 1D marcs model atmospheres (Gustafsson et al. 2008).Eleven elements are computed in non-LTE (Amarsi et al. 2020), the others in local thermodynamic equilibrium (LTE), and additional constraints on the surface gravities are incorporated through inferred distances by Bailer-Jones et al. (2021) and their limit on absolute magnitudes and luminosities.
In addition to the chemical information, we use the crossmatches of stars with value-added-catalog for stellar ages.The latter incorporates photometric and astrometric information from Gaia eDR3 (Lindegren et al. 2021) and 2MASS (Skrutskie et al. 2006) to estimate isochrone ages with the bstep code (Sharma et al. 2018).
We apply the following general quality and selection cuts to the data, firstly flag_repeat = 0 to get one measurement per star, and secondly flag_sp = 0, flag_fe_h = 0, flag_Mg_fe = 0, snr_c2_iraf > 25, and [Fe/H] > −2 to get the most reliable stellar parameters while maintaining a sufficient amount of stars.Whenever we use additional elemental abundances of an element X, we also enforce flag_X_fe = 0.
While we have considered enforcing bstep age uncertainties below 50%, we have decided against implementing a sharp limit, both to avoid complex influences3 in the selection function and because it does not significantly impact the stars with ages above 8 Gyr.
To allow better comparability with the simulations later on (see Sec. 3.1), we also limit our sample to stars with Galactic latitude || > 10 deg and distance   < 4.2 kpc (based on photogeometric distances from Bailer-Jones et al. 2021), the 95th distance percentile of GALAH DR3 (Buder et al. 2021).After applying all these cuts, our sample still consists of 284 923 stars, allowing a statistically meaningful comparison with theoretical predictions.

Theoretical predictions from a NIHAO Zoom-in simulation
The model perspective on accretion is provided through a cosmological zoom-in simulation of a Milky Way analogue (g8.26e11) from the Numerical Investigation of a Hundred Astronomical Objects (NI-HAO, Wang et al. 2015).The model galaxy has a total mass (gas, stars and dark matter) of 8.26 • 10 11 M ⊙ and contains 4 • 10 10 M ⊙ stellar mass with a stellar mass resolution of 1.06 • 10 5 M ⊙ (Buck et al. 2021) and was calculated as part of the NIHAO-UHD project (Buck et al. 2020).
Simulations were carried out with the smoothed particle hydrodynamics code Gasoline2 (Wadsley et al. 2017) using cosmological parameters from Planck Collaboration et al. (2014) with initial conditions and energetic feedback descriptions from the NIHAO project (Wang et al. 2015).Zoom-in simulations were then performed as described in detail by Buck et al. (2021) with star formation following Stinson et al. (2006) and energetic feedback following Stinson et al. (2013).
Because computational resources still limit the mass resolution of simulations, we are relying on tracer particles that represent simple stellar populations (SSPs) with the same age, overall metallicity and discrete initial mass function (IMF).Buck et al. (2021) have implemented the flexible chemical evolution code chempy (Rybizki et al. 2017) to calculate the chemical yields for the SSPs.In particular, we use the alternative (alt) setup of chempy that assumes a Chabrier (2003) IMF with high-mass slope of  IMF = −2.3 over a mass range of 0.1 − 100 M ⊙ for SSPs across a metallicity range of / ⊙ ∈ [10 −5 , 2].The code calculates the contribution from asymptotic giant branch (AGB) stars, CCSN across a mass range of 8 − 40 M ⊙ , and SNIa with a an exponential function with exponent −1.12, a delay time of 40 Myr, and a normalization of the SNI rate of -2.9.For each of these nucleosynthetic channels, yields from the following studies are used: Limongi & Chieffi (2018) for CCSN, Seitenzahl et al. (2013) for SNIa, and Karakas & Lugaro (2016) for AGB stars.
To achieve a roughly similar selection as the observational data (see Secs. 2.1 and 3.1), while maintaining enough star particles, we limit the simulated data set to star particles within a torus at 8.2 kpc galactocentric radius with a 4.2 kpc tube radius (corresponding to the 95th distance percentile of the GALAH sample) and neglect star particles within || < 10 deg or [Fe/H] < −2.We use the last simulation snapshot at a cosmic time of 13.8 Gyr to estimate the ages from the reported formation time.
The simulation traces the elemental abundances of H, He, C, N, O, Ne, Mg, Al, Si, P, S, V, Cr, Mn, Fe, Co, and Ba.Of these, ten correspond with the GALAH data and allow us to compare diagnostic plots possible with abundances of C, O, Mg, Al, Si, V, Cr, Mn, Co and Ba.To allow a better comparison of simulated to observed abundance patterns, we shift the abundance patterns [X/H] by the values listed in Table 1, such that the median abundance for [X/H] is Solar for the stars with ages of Solar age (4 − 5 Gyr) in the GALAH-like footprint of the Milky Way analogue.
The simulation records galactocentric Cartesian positions (,  , ) and velocities (  ,   ,   ) for each particle, which we while using the numpy arctan2 package.

DIAGNOSTIC PLOTS: OBSERVATION VS. SIMULATION
In Sec. 1, we have reviewed several observation-based plots that have been used to identify accreted stars in the Milky Way.In this section, we are now assessing if the simulated Milky Way analogue shows similar signatures when using these diagnostic plots.In Sec.1b when showing the extend of these stars with a convex hull outline.While this may seem surprising for observers at first sight, it has been seen previously also in the Auriga (Grand et al. 2020;Orkney et al. 2022), HESTIA (Khoperskov et al. 2023d, see their Fig.12), as well as the VINTERGATAN (Rey et al. 2023, see their Fig.A2) simulations, when tracing different merger subpopulations.
When applying a similar selection to NIHAO as would be expected from GALAH observations in the Solar vicinity (see Sec. 2.2), both the [Fe/H] vs. [Si/Fe] plane in Fig. 1d as well as the [Al/Fe] vs.
[Mg/Mn] plane in Fig. 1e cleared up significantly, with the accreted sequence becoming more prominent around [Al/Fe] of -0.3.
Intrigued by the difference in clarity of the chemical planes, we also show the two age-metallicity relations -or more accurately age-[Fe/H] planes -for the full Milky Way analogue in Fig. 1c and the The selection in panels a) and e) are drawn with a polygon with vertices at (-0.9,0.05),(-0.65,0.025),(-0.65,0.045),and (-0.9,0.07).A convex hull polygon of the selected stars is then drawn in other panels.
observational footprint in Fig. 1f.The clearly arising separations of sequences in both selections motivates a more dedicated analysis of this plane, as will follow in Sec. 4. The tight sequences, however, also raise the question, if our Milky Way analogue is indeed similar to our Milky Way, which will be discussed in Sec.5.3.Under this section's theme of identifying differences in the planes of full and limited sample selection, we notice a slightly larger scatter of abundances in Fig. 1c, that is, a differences of iron abundances outside the Solar vicinity, as well as the disappearance of a metal-poor ([Fe/H] ∼ −1.5) sequence below 10 Gyr when limiting the selection to the GALAH footprint.
Key Takeaway: Location matters, or better: a clever sample selection -even for simulations!

Abundance trends: [Fe/H] vs. [X/Fe]
In Fig. 2, we show the elemental abundances [X/Fe] as a function of iron abundance [Fe/H] of the 10 elements X that overlap between GALAH observations and NIHAO simulations.To first order we see an agreement of general trends for most elements, that is, elements that are produced by CCSN like O, Mg, and Si are showing a high plateau for low [Fe/H] and a decrease towards higher iron abundances -as expected by the onset of Fe-producing SNIa.The iron-peak elements Cr as well as Mn and Co show increasing abundances from low values at low [Fe/H] towards solar values at solar [Fe/H].Abundance estimates for V and Co do not agree, which is partially (but not fully) caused by unreliably high measurements of [V/Fe] and [Co/Fe] in GALAH (Buder et al. 2021).For Ba, we note tightly increasing enhancements towards Solar [Fe/H], starting from initially low values around [Ba/Fe] ∼ −0.5.Furthermore, we notice a clear second (accreted) sequences below [Fe/H] < −0.5 for C as well as Ba.While not shown in this manuscript, we note that when plotting the full simulated galaxy in Fig. 2, no significant difference in trends can be identified, but a larger spread of abundances -most notably towards the lowest metallicities.
When focusing on the actual absolute abundances (note the difference in y-axis), we note significant differences between the chemical evolution sequences.To understand this difference, we first need to remind ourselves that the observational data is limited by the ability to measure abundances (note the missing abundances of [C/Fe] for low metallicities due to detection limits) and includes measurement noise, while the theoretical data does not necessarily have to reproduce the formation scenario of the Milky Way and chemical enrichment tracing is limited to our best current understanding of nucleosynthesis, see e.g. the discussion of nucleosynthetic yields in Buck et al. (2021).
That said, without adjusting for abundance-offsets, we note that none of the abundance sequences in the Solar metallicity regime of [Fe/H] ∼ 0 intercept with the expected Solar value.While C, O, Si, Cr, Mn, and Ba are significantly overproduced between 0.05 (for Cr) and 0.5 dex (for Ba), we notice a significant under-abundance on the −0.3.. − 0.1 dex level for Mg, Al, V, and Co.Nevertheless, for the remainder of this manuscript this poses no severe limitation as we are mainly interested in a differential signal in order to separate in-situ and accreted stars.
Key Takeaway: While simple abundance-metallicity plots do not easily reveal accreted structures for most elements, both C and Ba abundances show the potential to identify different abundance sequences -if the nucleosynthetic processes of the simulation are accurate.

Abundance-abundance plots: [Al/Fe] vs. [Mg/Mn]
We now turn to abundance-abundance plots of [Al/Fe] vs. [Mg/Mn] that were already mentioned in Sec.3.1.These plots combine the enrichment channels of CCSN vs. SNIa, that is [Mg/Mn], and the enrichment of Al with more complex metallicity-dependent enrichment pathways (Hawkins et al. 2015;Das et al. 2020;Kobayashi et al. 2020).The enhancements of these elements should differ significantly when comparing in-situ stars born in the early Milky Way with the less massive system that it accreted because of the expected lower star formation intensity (and thus iron abundance) and occurrence of CCSN in the latter.
When using the same axis limits for simulation and observation (Figs.3a and c), we notice the strong offset and significantly smaller spread of abundances in the simulation, in agreement with our pre-  vious findings for abundances in Sec.3.2.When zooming into the relevant [Al/Fe] vs. [Mg/Mn] regime for the simulations, however, we notice a better qualitative agreement, in that we see two significant overdensities for simulated and observed galaxy.For the observations, we identify the lower and upper overdensity at [Mg/Mn] of 0.0 and 0.4.These have previously been confirmed as low-and highα-sequence (or thin and thick disk), respectively (Hawkins et al. 2015;Buder et al. 2022).These sequences show similar Solar or even enhanced [Na/Fe] and [Al/Fe].In contrast to this, the major overdensity in the simulation is found around (-0.25,-0.5)with tails reaching higher [Mg/Mn].The minor overdensity is found towards higher [Mg/Mn] as well as lower [Al/Fe] -coinciding with the expected region of accreted stars (Horta et al. 2021) which are less numerous in the observational data of GALAH DR3.
To achieve a more realistic comparison, we add noise on the order of 0.01 and 0.04 dex to the [Al/Fe] vs. [Mg/Mn] plane of the NIHAO simulation in Fig. 4. It has to be noted though, that this noise has to be compared to the separation of accreted and in-situ sequence, which is on the order of 0.07 dex for [Al/Fe] in the simulation (see black dashed lines in Fig. 4), compared to 0.18 dex in the Milky Way (Buder et al. 2022, see their Tab. 4).This means that a noise of 0.03 dex (Δ ∼ 2.3) in NIHAO has a similar effect as the noise of 0.08 dex in GALAH DR3.When comparing the different panels, we see that the initially very clear separation of the accreted and in-situ sequence quickly washes out to the point of current observational noise (between Figs.4c and 4d), where more sophisticated methods are needed to assign particles to either sequence, as done in the past by observers (Das et al. 2020;Buder et al. 2022).
Key Takeaway: The abundance-abundance plots of [Al/Fe] vs. [Mg/Mn] show qualitative agreement with the observational data, but with significant abundance offsets and a more pronounced overdensity of accreted stars than the observations in the absence of noise.This could indicate a more massive merger progenitor than the Milky Way's GSE or hint at our inability to actually measure Mg, Mn, and Al or Na abundances of accreted stars as well as reflect theoretical uncertainties in the yields of Mg, Mn, Al and Na.When adding realistic uncertainties to this plane, the separation washes out to the point where an assignment to either sequence needs more advanced methods.

A NEW ANGLE: AGE-ABUNDANCE-DISTRIBUTIONS
Fortunately, chemistry is not the only conserved property of stars.We therefore turn to stellar ages and assess age-abundance trends in the data, starting from the age-[Fe/H] relationship that we already showed in Figs.1c and f.This relation has proven to be a powerful tool not only for resolving processes with observed data of our Milky Way (e.g.Twarog 1980;Edvardsson et al. 1993;Nordström et al. 2004;Casagrande et al. 2011;Feuillet et al. 2019;Xiang & Rix 2022), but also other galaxies -for example to identify accreted stars (e.g.Pinna et al. 2019a;Martig et al. 2021).However, observations still suffer from larger uncertainties on the order of 15 − 50% for most stars, so examining simulation plots removes this barrier and allow us to determine chemical enrichment as a function of age.Additionally, turning towards stellar ages in the simulations focuses on the fundamental variable age and removes uncertainties associated with uncertain yield sets.

Clearly separated sequences in the age-[Fe/H] relation of the Milky Way analogue
Based on Fig. 1f, we identify two definite and one tentative sequence in the age-metallicity relation with the following selections, selecting 75 440 particles for Sequence 1, 13 398 particles for Sequence 2, and 133 particles for Sequence 3 of the 88 980 particles in the footprint: We discuss these numbers in comparison with the Milky Way later in Sec.5.3.In this section, we are primarily interested in tracing the three sequences across their chronochemodynamic properties, as done with separate rows in Fig. 5, starting with the age-[Fe/H] relation.
For Sequence 1, the most metal-rich sequence with 85% of particles in the footprint, we see the stars being distributed closest to the galactic plane (| | < 2 kpc) in Fig. 5b.Their kinematics in Fig. 5c) likens a Toomre diagram of the Milky Way (Helmi et al. 2018, see their Fig.1a) but with a smoother transition of thin and thick disk.The smooth transition of thin and thick disk in kinematic space could lead to the question if the Milky Way analogue truly has a thick disk (see Bovy et al. 2012, who raised the same question for the actual Milky Way).This impression is confirmed in chemical space of [Fe/H] vs.
[Mg/Mn] plane.We refrain from further dissecting this population in old and young sequence to establish a thin and thick disk as this is not the focus of this work.
Sequence 2 with 15% of particles in the footprint shows the behaviour of kinematic halo stars with a spatial distribution not restricted to the galactic plane (Fig. 5g).While Buder et al. (2022) found   = 20 +100 −70 kpc km s −1 when selecting the major accretion remnant in the Milky Way via chemical compositions, we find significantly prograde orbits with   = 150 +100 −100 kpc km s −1 (Fig. 5h), suggesting a net rotation almost as large as the rotation of sequence 1 (  = 180 +60 −80 kpc km s −1 )4 .Once again, we notice the offset towards lower abundance enhancement of this sequence in the [Fe/H] vs. [Si/Fe] plane as well as the [Al/Fe] vs. [Mg/Mn] plane with respect to sequence 1 (Figs.5i and j, respectively).
Finally, we notice that the smallest Sequence 3 (only 0.1% of the stars in the footprint) significantly overlaps with Sequence 2 in all of the chemodynamic properties.In addition to the very clear separation in age-abundance space (Fig. 5k), the only way to identify particles of this sequence would be through their relative overdensities in the abundance-abundance planes (Figs.5n and o).
Key Takeaway: We find a smooth transition of thin and thick disk (identified as Sequence 1) of the Milky Way analogue in Figs.5a-e.Sequence 2 clearly presents itself as an accreted sequence in any of the plots of Figs.5f-j.In particular, the dominance of the thin/thick disk makes it hard to identify halo stars in spatial and kinematic space, wheres age-abundance and abundance-abundance plots show less contamination.Finally, we note that the smallest sequence (with lowest [Fe/H] values) is overlapping with the major accreted sequence spatially, kinematically, and chemically, and can only be separated as clump in the age-abundance (Fig. 5k) or abundance-abundance space (Figs.5n and o).Having established the age-abundance planes as the best diagnostic plane (in the absence of measurement uncertainties), we are now concerned with the information held in age-abundance relations beyond the age-[Fe/H] relation.

Age-abundance relations beyond age-[Fe/H]
Given the various abundance distributions in the [Fe/H] vs. [X/Fe] space from Fig. 2, we are now interested in the trends for different elements with age.Previous research has suggested that not necessarily [X/Fe], but [X/H] may be a better tracer of accreted stars (e.g.Fuhrmann et al. 2017;Feuillet et al. 2021, for [Mg/H]).Furthermore, we are inspired by the work by Weinberg et al. (2019Weinberg et al. ( , 2021) ) as well as Griffith et al. (2019Griffith et al. ( , 2022)), who modeled nucleosynthetic contributions by SNIa and CCSN with amplitudes such as  cc = 10 [Mg/H] ( Weinberg et al. 2019).To avoid confusion with both the abundance notation (X) and their notation  cc / Ia , we want to explicitly lay out our notation of rewriting the element abundances in a linear way as  X/H , so that they are linearly proportional to the ratio of number densities  X / H :

Age-[X/Fe] or age-[X/H] or age-𝑁 X/H ?
In Fig. 6 we show the age-abundance distributions for six of the elements overlapping between NIHAO and GALAH, namely C, O, Mg, Al, Mn, and Ba.The other abundance trends are appended in Fig. A2 as their trends liken those already shown in Fig. 6 and they thus do not add new insights.
First, we discuss the top panels of Fig. 6 showing the age-[X/Fe] trends.As pointed out in our description of Fig. 2, the NIHAO simulation show a rather flat, enhanced [C/Fe] abundance for accreted stars, that is, [C/Fe] ∼ 0.11 compared to [C/Fe] < 0.07 at around 10 Gyr.For O in Fig. 6b, we notice a decrease of high [O/Fe] abundance with decreasing stellar age and an emerging bifurcation with higher [O/Fe] ∼ 0.1 for accreted stars.For Mg, the in-situ and ac-creted components overlap in the age-[Mg/Fe] plane (Fig. 6c).For both [Al/Fe] and [Mn/Fe], we see 0.05 dex lower abundance for accreted stars than for in-situ stars, with a clearer separation of sequences for [Al/Fe] at oldest ages.Finally, the [Ba/Fe] plane shows a lower initial abundance of [Ba/Fe] at oldest ages, which is expected -at least qualitatively -due to the delay of Ba-production by AGB stars (Karakas & Lugaro 2016).We note here that the production of Ba relative to Fe for accreted stars actually surpasses those of in-situ ones, as can be seen from the crossing of Sequence 1 and 2 in Fig. 6f at around 9 Gyr.At face value, these differences would exclude both the age-[Mg/Fe] and age-[Ba/Fe] plane as useful diagnostics, whereas all other elements show offsets in the age-[X/Fe] plane for at least the youngest (8 − 11 Gyr) accreted stars.Moving towards the age-[X/H] plane in the middle row of Fig. 6, we notice a very similar qualitative behaviour of all elements.However, we note that the separation of accreted and in-situ sequences in this plane is roughly five times larger, that is, around 0.5 dex, than in the age-[X/Fe] plane.Finally, we consider the age- X/H plane in the bottom rows of Fig. 6, where we again see a qualitatively similar trend.When focusing on the accreted sequence, we note a roughly linear enhancement of  X/H for all elements except for Ba, which shows an increasingly enhanced curve towards younger ages, in line with an expected increase of Ba-production by AGB stars (Karakas & Lugaro 2016).The difference between sequences starts very small at oldest ages, but then increases to Δ X/H ∼ 0.7 around 8 Gyr.While it is not the focus of this particular work, we want to note the increase in scatter of the in-situ age-abundance sequence, most notably in  X/H , for in-situ stars younger than the oldest accreted stars.This pattern has already been recognised in previous work on NIHAO-UHD simulations (e.g.Lu et al. 2022) and was connected with the time of disc formation for several Milky Way analogues.While Lu et al. (2022) identified the strong correlation with disc formation time (see their red lines in their Fig.1), they also noted a correlation with major mergers (Lu et al. 2022, see their orange lines in Fig. 1).The analysis of different simulation snapshots is beyond the scope of this work.But for completeness we want to note the strikingly coinciding times of the onset of scatter in the in-situ sequence and end of star formation in the accreted sequence for the NIHAO-UHD Milky Way analogue of our study, g8.26e11 which signifies the point in time of accreting the dwarf galaxy progenitor onto the Milky Way.
Returning to the main focus of this work, that is, to better identify accreted stars in observations, we have created Fig. 7 with observational data from the GALAH survey (Sec.2.1).While we see similarities between observations and simulated age-abundance trends for several elements, other elements show different behaviour.In particular, the measurements of C in GALAH are less numerous for older stars.This is expected based on the difficulty of measuring C from the atomic CI line.We will therefore refrain from a comparison for this element.For O, Mg, and Al, we see an increase of [X/Fe] towards old ages, For Al (Fig. 7d), this is contrary to the trend of the simulations (Fig. 6d).The plots with respect to [X/H] are rather flat and inconclusive.For Mn, we see a decrease of both [Mn/Fe] and [Mn/H] towards old ages, in agreement with the simulation.For Ba, we see a decrease of both [Ba/Fe] and [Ba/H] towards old ages -in agreement with the simulations, but even more interestingly the hint of an overdensity of high [Ba/Fe] ∼ 0.2 and low [Ba/H] ∼ −0.6 around 9 Gyr.This unique behavior is also seen in the simulations -and also only for Mn.When looking at the age- X/H plots in the lower row of Fig. 7, we identify the aforementioned overdensities again for  Ba/H at 9 Gyr, but also  Mn/H (Fig. 7r), whereas the other elements C, O, Mg, and Al do not reach the regime of  X/H < 0.4 in larger numbers.
Key Takeaway: Across the different age bins, we find Gaussian distributions of abundances in the [X/Fe] and [X/H] space.At face value, the abundance vs age plots are separated most clearly when the abundance are noted relative to H rather than Fe.Establishing the relationships relative to the number density  / should also be explored in future work as they bear the potential to allow the tracing of linear and non-linear contributions through chemical yields.

Quantifying the difference of age-[X/H] sequences
In this section, we are concerned with the significance of separation between the in-situ and accreted sequences in the simulations -for now in the idealised case without observational noise (we will study the latter in Sec.4.2.3).
For this purpose, we trace the distribution of abundances [X/H] for the three sequences defined by Eqs.4-6 in histograms across different age bins.In Fig. 8, we show the histograms for the different sequences in age bins of 0.5 Gyr above 8.5 Gyr, the youngest age of stars in Sequence 2. The panels of each figure include further useful information.Following the color-scheme from Fig. 5, they state the relative contribution of each sequence to the age bin in percent, followed by the median abundance as well as the difference to the 16 th and 84 th percentile.We note firstly that the contribution of accreted stars is increasing from 8% for 8.5 − 9.0 Gyr to 78% for 12.5 − 13.0 Gyr.In addition to the decrease of overall [Mn/H] towards older ages, we note that the scatter of both distributions   increases from the youngest stars (∼ 0.02 − 0.05 dex) to the oldest (∼ 0.1 − 0.4 dex), while the distributions remain in a Gaussian shape for the most part, with larger tails towards lower [Mn/H] only noticeable for the oldest stars (> 12.5 Gyr).
To determine the significance of separation, we are inspired by the previous study by Lindegren & Feltzing (2013), who computed a separation significance  based on the differences of mean abundances weighted by the combined standard deviation of abundances.Having convinced ourselves that the distributions are Gaussian in the previous section, we now expand their formulation by adding observational Gaussian noise as another additive to the weights and thus arrive at the formulation Assuming no observational noise, we calculate  12 based only on the means and standard deviations for the age bins and elements [X/H] listed in Tab.A1 for the sequences defined in Eqs.4-Eqs.6.For completeness, we also list the properties of Sequence 3, when more than 100 particles are present.
Key Takeaway: The intrinsic scatter in [X/H] for both in-situ and accreted sequences is larger (∼ 0.1 − 0.4 dex) for older stars (> 12 Gyr) as a result of larger change in this logarithmic scale when compared to stars born at the younger age of 8.5 − 11.0 Gyr (∼ 0.02 − 0.05 dex).This would make the younger accreted stars a better candidate to trace the chemical enrichment history and infer progenitor masses than old stars.While most elements, with the exception of the less suitable Ba, showed a similarly large separation significance, the odd-Z elements P and Al would be slightly more preferable than the iron-peak and alpha-process elements, which all outperform the lightest simulated elements C, N, O, and Ne.

A more realistic picture: The influence of observational noise on age-[X/H] sequence
Having determined the most suitable plots for identifying accreted populations in simulations, we now consider this problem from a more realistic perspective and quantitatively determine how much noise we could allow for both abundance and age estimates to still be able to separate the sequences.
To get an intuition for the influence of measurement noise on the separation significance of abundances, we test different ranges of abundance uncertainties of  obs as input for Eq. 8.The resulting values are visualised in Fig. 9 and representative values listed in Tab.A2, which also includes median uncertainties of the measurements for GALAH DR3 (Buder et al. 2021), ranging from  GALAH = 0.06 for [Si/H] to  GALAH = 0.13 for [O/H] and [V/H].
We expect the sequences in age-[Fe/H] will be significantly more difficult to distinguish in older stars with corresponding larger uncertainties.To get a better intuition, at which age uncertainties we gain or loose more diagnostic power, we are performing a simple experiment.Assume that we know that our best chances to distinguish the sequences is to look at the youngest accreted stars, that is, those within our data set closest to 8.5 Gyr.Using the smallest reasonable age bin (8.5−9 Gyr), we estimate a reasonable separation of Sequence 2 to 1 by estimating the mean and standard deviation of the sequence and adding a reasonable observational noise of  obs = 0.05 dex to it.This allows us to discern Sequence 1 from 2 with a rather simple cut of stars being above or below  2 + 3 √︃  2 2 +  2 obs , as shown in Fig. 10a for [Fe/H] ≶ −0.49.Using this abundance limit as the boundary, we can then test the distributions of Sequence 1 and 2 for increasingly large age bins of 0.5..(0.5)..4.0 Gyr, as visualised with the selection boxes in Fig. 10a.Performing this experiment for the various elements, we reach the distributions of  12 in Figs.10a and 10b with increasing age bins size (which we equate with twice the age uncertainty).
Because we are working within the age range of 8.5−12.5 Gyr, this means that an age bin of 2 Gyr would be equivalent to an age uncertainty of just over 20%.We can see from the plot that we require an uncertainty of less than 20% to achieve the highest separation significance -value.While being an insightful gedankenexperiment, our tests on the influence of abundance and age uncertainties, however, neglect an important fact.That is, that both uncertainties are present at the same time in observations in complex ways and a selection like ours with prior knowledge of the sequence positions is not possible for the Milky Way.
Guided by our previous results on abundance and age uncertainties, and with the actual uncertainties of GALAH DR3 in mind, we have added both noises to the properties of the simulated galaxy and replicated the age-[Fe/H] in Fig. 11.This idea was also neatly implemented by Renaud et al. (2021a, see their Fig.12) to create more realistic noisy data for the VINTERGATAN simulations.Our figure shows the extremes of no noise on the far left (Fig. 11a) and the actual observed Milky Way age-[Fe/H] relation in the Solar neighbourhood on the far right (Fig. 11e with an average age uncertainty of ∼ 25% and [Fe/H] uncertainty of 0.075 dex).In between, we show different realisations of added Gaussian noise to the simulations.The chosen noise levels mimic those reported by pioneering studies, such as the one by Xiang & Rix (2022) with a reported median uncertainty of 7.5% for age and 0.05 dex for [Fe/H] (see Fig. 11b), the series of studies by Nissen & Schuster (2010) and Schuster et al. (2012) reporting uncertainties of 0.03 − 0.04 dex for [Fe/H] and a median of 15% for age (Fig. 11d), as well as the already mentioned GALAH DR3 age-[Fe/H] relation (compare Figs. 11d and 11e) with uncertainties of around 25% for age and 0.075 dex for [Fe/H].
We want to stress here, that we assume that the uncertainties of metal-rich and metal-poor as well as young and old stars are relatively the same.In reality, this is not the case, because the measurements of old and metal-poor stars are less certain (Frebel & Norris 2015) when compared to stars with similar ages and metal-content than our Sun.The estimation of uncertainties -most importantly accuracy -is complicated.Most studies and surveys thus only report precision uncertainties, which neglects uncertainties of stellar evolutionary models (Soderblom 2010) just as much as uncertainties of abundance inference from stellar spectra, for example with approximations made in the modelling of atmospheres with uncertain line transitions (Mihalas & Athay 1973;Asplund 2005), among other factors (Jofré et al. 2017(Jofré et al. , 2018;;Nissen & Gustafsson 2018).
Setting these concerns aside when looking at Fig. 11, we can easily see that a larger noise in either direction blur out the clearly separated sequences of the noiseless simulation.While we can still see the separation of sequences in Figs.11b and 11c for the youngest accreted stars, the sequences overlap too much in Fig. 11d to be able to tell them apart by eye.We would therefore not expect to be able to tell apart a similarly separated sequence in the observed age-[Fe/H] of GALAH DR3 (Fig. 11e) in this particular plane.It is therefore not surprising, that the now established substructures of the stellar   4) and later by Xiang & Rix (2022) who each looked at main-sequence turnoff and subgiant stars, for which observable properties change more significantly as a function of stellar age.In particular, Xiang & Rix (2022) were able to resolve a "Z-shaped structure" in the age-[Fe/H] relation of stars with lower angular momentum (see their Fig. 2).We note that the third sequence (bottom of Fig. 5) would not be detectable in this plane even with the currently lowest uncertainties.
For completeness, we also want to note that the simulation at similar noise levels to GALAH DR3 shows a much tighter disk sequence than the actual Milky Way (compare the top of Figs.11d and 11e).This abundance spread would only be replicated if we were to add a noise of  [Fe/H] > 0.15 dex to the simulation.This either suggests that the well-documented radial migration (e.g.Haywood 2008;Frankel et al. 2018) of more metal-rich stars from the inner Galaxy and more metal-poor stars from the outer Galaxy to the Solar vicinity is not happening at the same strength in the Milky Way analogue or that turbulent metal mixing is slightly underestimated in the NIHAO simulations (see also discussion in Buck et al. 2021).
Key Takeaway: Testing the uncertainties in the age-abundance planes independently and in concert, we find critical uncertainties of  [/ ] = 0.15 dex and  age = 20% above which we loose the ability to tell apart accreted and in-situ sequences with certainty.Adding noise at the level of current pioneering studies to the simulations confirms the qualitative impression of a "Z-shaped structure" (Xiang & Rix 2022) in the noisy age-[Fe/H] plane for two underlying sequences (Fig. 11b).

DISCUSSION
The initial motivation of this research was to create a variety of abundance planes and compare them using both the simulated and observed data to determine the best abundance-abundance plots for identifying substructure that could be representative of an accreted population of stars.In this section, we therefore discuss not only the question which information reveals accreted structures best (Sec.5.1), but also take a closer look at the age-[Fe/H] relation around the time of the last merger and discuss implications of the found patterns (Sec.5.2).At the end of our discussion, we touch on how comparable the simulation at hand is to the Milky Wayand what that means for the findings of our comparison; also in the context of galaxy evolution (Sec.5.3).

Which information reveals accreted structures best?
In an evolving galaxy, the kinematic and dynamic information of stars that were born together is expected to change.This memory may not be completely lost throughout a galaxy's vivid evolution.This is exemplified by the many recovered structures (Deason & Belokurov 2024) and even metallicity gradients (Khoperskov et al. 2023a) in velocity and energy-angular momentum space.Nonetheless, orbits are certainly expected to be blurred throughout evolution.Even within our simulated data, we see significant overlap in the kinematic characteristics of the three separate sequences (Fig. 5).Depending on the size and time of mergers there could be significantly more overlap between accreted and in-situ populations.
Chemical compositions and ages -although harder to infer -are expected to be more conserved in stellar atmospheres than dynamic properties, a paradigm of the idea of chemical tagging (Freeman & Bland-Hawthorn 2002).Our analyses, in agreement with previous analyses (e.g.Khoperskov et al. 2023d;Rey et al. 2023), however found that while stars that were born together may share similar chemistry, this chemistry does not have to be significantly different from other birth places for the light elements inspected in this study.This problem brings the study by Ness et al. (2018) to mind, who found it difficult to tell apart stars of the stellar disk with light element abundances for 0.3% of their inspected field star pairs.Extending this study to heavy elements with GALAH, (Manea et al. 2023) found that heavy elements help to decrease such doppelgänger rates by an additional factor of 6.The simulation of this study does only follow up Ba with the currently limited understanding of neutron-capture enhancement.Although we have highlighted the interesting observational pattern of age-[Ba/H] and age- Ba/H earlier (see Sec. 4.2.1 and the discussion of Figs.7f, 7l, and 7s), this avenue needs to be explored in future studies -both observationally and through chemical evolution modelling.
In some cases the chemical information of some structures (e.g.Sequence 3) is quite compact in abundance-abundance space (see Fig. 5n for [Fe/H] vs. [Si/Fe]).A careful selection and combination may still provide a useful avenue to trace back such substructure that has lost its kinematic memory, even though not all chemical elements may provide useful insight (see e.g. for Al in Fig. 5o).
Again, it is important to note that the simulated galaxy may have different characteristics to our own Galaxy.As discussed, there are two clear streams that have been attributed as accreted and in situ, based on the assumption that the accreted stars would have been born in a smaller galaxy with less massive stars and therefore lower element abundances.
While the jury on the best information for identifying is still out, this study yet again underlines that the survey of the whole galaxy is important to identify smaller accretion events.While for example a non-negligible amount of stars of the accreted Sequence 2 can be found in the observable footprint (see Fig. 12a), we find stars of Sequence 3 only far outside of the GALAH-like footprint at galactic heights of || > 10 kpc (see Fig. 12b).This area is being probed by the H3 Survey (Conroy et al. 2019) in the actual Milky Way and has indeed found a significant amount of substructure, with more findings expected by the upcoming 4MOST halo surveys (Christlieb et al. 2019;Helmi et al. 2019) as well as the 4MOST dwarf galaxies and streams survey 4DWARFS (Skúladóttir et al. 2023).

A closer look at age-[Fe/H] relation around the last major merger
While stellar mass estimates can help us to quantify the mass ratios of previous mergers, we ultimately want to probe not just the presence and correlation but the causality of major mergers and chemical thin/thick disk bimodality, among other quantities.To this end, we are thus focusing on the age-[Fe/H] relation around the time of the last major merger in this section.While this is important both in terms of stellar particles and gas particles (see e.g. the study by Buck et al. 2023), this study focuses on the stellar particles as the surviving and observable tracers throughout galactic evolution.In Fig. 13, we are thus zooming into the age-[Fe/H] relation around 8.6±0.4Gyr, the rough 5 time of the last merger.Around this time, we notice not only the stop of star formation in the accreted sequence (see 5 Quantifying the beginning and end of a merger process with all its intricate dynamical and chemical interactions is difficult.However, our rough estimate of 8.6 ± 0.4 Gyr is almost twice as precise as current best age precision uncertainties. light blue region in Fig. 13a), but also the appearance of several particles with iron abundance between the in-situ and accreted sequence (see pink region in Fig. 13a).We thus trace several properties of the in-situ (dark and light orange) as well as accreted sequence (dark and light blue) before and throughout the merger, respectively.In particular, we notice a continuing increase of [Fe/H] between 9 − 9.5 Gyr to 8 − 8.5 Gyr for both sequences.For the in-situ sequence, we notice an increase from −0.165 ± 0.038 to −0.113 ± 0.056.In particular, we find an increase of the [Fe/H] spread between 9 − 9.5 Gyr to 8 − 8.5 Gyr from 0.038 to 0.056, that is an almost 50% increase.When calculating the same quantities for the whole galaxy (Fig. A3), we find a smaller increase of [Fe/H] spread from 0.065 to 0.076, that is only 17%.For the accreted sequence, we find an increase from −0.682 ± 0.029 to −0.634 ± 0.031 with consistent scatter.We further find a decrease in the ratios of [Mg/Fe] and [Mg/Mn] in Figs.13b  and 13c, whereas the ratio of [Al/Fe] stays rather constant for both sequences (Fig. 13c).
Most interesting to us seem the pink particles of Fig. 13 with [Fe/H] between the in-situ and accreted sequences.These arise throughout the merger with ages of 8.6 ± 0.4 Gyr and are also exhibiting intermediate abundances in [Mg/Fe], [Al/Fe], and [Mg/Mn] in Figs.13b and 13c.Such stars have previously been identified as stars that formed further outside of the Solar vicinity with larger angular momenta (e.g.Haywood 2008;Haywood et al. 2013;Buder et al. 2019).We too find them to have larger angular momenta of   = 1340 680 −490 kpc km s −1 when compared to   = 1290 600 −480 kpc km s −1 for stars of the in-situ disk at the same time and even higher than the accreted stars with same ages (  = 1140 700 −600 kpc km s −1 ).The relative contribution of stars with these properties is, however, relative little with a stellar mass of only 0.1 × 10 8 M ⊙ compared to the amount of stars born in the dark and bright orange in-situ regions (6.7 × 10 8 M ⊙ and 6.1 × 10 8 M ⊙ ) and even the accreted stars with their decreasing star formation between dark and bright blue regions (1.2 × 10 8 M ⊙ and 0.4 × 10 8 M ⊙ ).Even when extending the analysis towards younger in-situ stars with lower [Fe/H] (see Fig. A4) we continue to find relative low numbers of these stars, whereas they seem to be more numerous in the Milky Way's Solar vicinity (Haywood et al. 2013;Hayden et al. 2015;Buder et al. 2019), thus raising important question of the similarity of Milky Way and simulated analogues which we discuss in the next subsection.

How comparable is the NIHAO Milky Way analogue to the
Galaxy?
Throughout this study we have tried to quantify the usefulness of specific diagnostic plots based on noiseless simulations of a Milky Way analogue.In this section of the discussion, we are now concerned with the question, how similar the Milky Way analogue is to the Milky Way and to which extend our conclusions can thus be applied to the Milky Way.

Chemical comparability:
As we acknowledged throughout this study (see e.g.Sec.2.2), the chemical evolution of Milky Way analogue simulations is neither fully accurate nor consistent with the Milky Way.While the basic properties of the galaxy at hand agree with the Milky Way within the order of magnitude, we acknowledge that a different merger history may lead to different chemical compositions.More of concern for our study is the inaccuracy and incompleteness of chemical enrichment parameters, such as chemical yields.While we can overcome potential relative inaccuracies in element production by up-or down-scaling the total number densities of elements (which in a logarithmic abundance notation leads to simple  shifts, as tabulated in Tab.1), we cannot correct missing channels -for example for Ba or other heavy elements.This certainly means that the chemical evolution and predicted sequence separations have to be taken with a grain of salt and can only be used as indicators.
Merger mass ratio: For the Milky Way, different estimates of the last major merger have been estimated through different comparisons.
When tracing the stellar masses of the in-situ and accreted sequences (Fig. 14), we find a similar stellar mass ratio of 1:3.0 around 8.2−8.6 Gyr based on stellar masses of 0.8×10 9 M ⊙ and 2.4×10 9 M ⊙ for stars older than 8.6 Gyr.
The (missing?) separation of thin and thick disk: Contrary to the more distinct sequences that were found by Nissen et al. (2020); Sahlholdt et al. (2022); Xiang & Rix (2022) for thin and thick disk stars in the age-[Fe/H] relation, we do not see a significant change along the age-[Fe/H] relation of in-situ component (as in thin vs. thick disk).Similarly, the [Fe/H] vs. [Mg/Fe] sequence does not present itself as the clear bimodal one seen in the Milky Way (Hayden et al. 2015) and other NIHAO simulations (Buck 2020;Buck et al. 2023).This has several reasons: a) the Milky Way analogue discussed in this study uses a different yield set (alt as shown in Buck et al. 2021) but traces many more elements.Comparing to Fig. 9 of Buck et al. (2021) we can see that this yield set produces a tighter sequence not quite resembling the observed one.b) the particle mass of the simulation is 10 5 M ⊙ for star particles which for numerical reasons diminishes metal mixing and a more clearer separation of the sequences.Additionally this means we are not resolving all merging dwarf galaxies but only the most massive ones.Nevertheless, the approach taken in this study by only looking at differential analyses of abundance patterns are robust findings.The presence of a thick disk may not necessarily be seen as a clear separation in [Fe/H] vs. alpha-enhancement, despite manifesting itself in the kinematics (see e.g.Minchev et al. 2013, their Fig.12).
This Milky Way analogue did only undergo one significant merger around 8.5 Gyr ago.We find an absence of a clearly chemically distinct disk with a more or less significant interruption in star formation and change in later chemical evolution of disk stars.A similar finding was made by Minchev et al. (2013) where in the absence of early-on massive mergers, their thick disk was kinematically less hot with vertical velocity dispersions smaller by a factor of ∼ 2 when compared to the Milky Way.Extending their conclusion into the chemical dimension, our analysis of only one(!) Milky Way analogue could lead to a similar conclusion as theirs, that is, that "the Milky Way thick disk is unlikely to have been formed through a quiescent disk evolution" (Minchev et al. 2013).
The (missing?) spread in the age-[Fe/H] relation of the young disk: In Sec.5.2, we pointed out the low spread of [Fe/H] among young disk stars.From 9.0-9.5, 8.5-9.0, and 8.0 − 8.5 Gyr, we see a change of [Fe/H] from −0.166 ± 0.037 to −0.115 ± 0.037 to −0.109 ± 0.050 in the simulation.This compares to a change of −0.222 ± 0.255 to −0.203 ± 0.258 to −0.186 ± 0.256 for GALAH DR3, respectively 6 .This significantly lower spread of the simulation could either specifically point towards a different formation history of disk stars in our simulation because of a much lower radial migration than in the actual Milky Way (Minchev & Famaey 2010;Loebman et al. 2011;Frankel et al. 2018) or reflect a suppressed metal mixing in the simulations.

Bottom line:
We acknowledge that the discussion of this particular section is incomplete and more comparisons could be done.However, we believe we have established that in terms of accretion, the simulation at hand can provide useful insight and we therefore conclude that further discussion about the comparability of simulation and observation (such as gas masses or the presence and influence of a bar or spiral arms) are important next steps but beyond the scope of this work.Such research should be performed also with other cosmological simulations, preferably with varied chemistry and formation histories.This would enable researchers to observe exactly how different variables in the simulation affect outcomes.Such comparisons would not only be insightful for the Milky Way, but also observations of other galaxies (compare for example figures by Pinna et al. 2019b,a;Martig et al. 2021, to our Fig. 12) and the ongoing observational projects such as the large MUSE program GECKOS (van de Sande et al. 2023).

CONCLUSIONS
Both the research field of Galactic archaeology in particular and galaxy evolution as a whole aim to unravel how galaxies form and evolve across cosmic time.While the observation of galaxies across different redshifts can provide us with insights into the pathways of galaxy evolution, the Milky Way provides us with the unique resolved information at redshift zero.To link all of these different snapshots, we can use the numerous scenarios provided by simulations, which now also allow us to trace the chemical evolution across different elements based on our current best understanding of nucleosynthesis.Mergers of galaxies and accretion of smaller systems are common 6 Here we report median values and half the difference between 16 th and 84 th percentiles.Assuming a Gaussian distribution, we would measure larger scatters with a qualitatively similar behaviour for the overall [Fe/H] values of −0.165 ± 0.038, −0.114 ± 0.043, and −0.113 ± 0.056 for the simulation and −0.251 ± 0.292, −0.225 ± 0.281, and −0.205 ± 0.271 for GALAH DR3.
events in a Universe that forms hierarchically and we see evidence of such mergers especially at high redshift  > 2 (or more than 10 Gyr ago).
An outstanding question is, if mergers play a minor role or major role in shaping galaxies.Do they for example only drive the scatter in the relation of age (or redshift) with velocity dispersion and metallicity, or do they actually set these relations and influence the shape of galaxies through the triggering of disk and bar formation?
To get a quantitative handle on answering these questions, we need to trace back the merger history of the Milky Way as the one galaxy for which we can get a resolved understanding of evolution across time and space.This paper thus aims to find new insights from the cosmological NIHAO simulations on how we can identify the surviving stars of accreted galaxies.
The major take-away points from this undertaking were: • A first assessment of the common diagnostic plots of observations manifested the old credo that the absence of clear features or subpopulations in the data does not mean that they are not presentlocation and selection do matter (Sec.3.1).While the assessment of chemical abundance planes of the Milky Way analogue for the whole galaxy is actually overshadowing the minor contribution of accretion, we confirmed the potential of age-abundance and abundanceabundance plots for smaller regions, both in similar footprint as current Milky Way observations (Fig. 1) as well as extragalactic observations (Fig. 12).
• In particular, when looking at abundance with respect to iron abundance (Sec.3.2 and Fig. 2) as well as [Al/Fe] vs. [Mg/Mn] (Sec.3.3 and Fig. 3) we confirmed the strong potential of chemistry to identify major substructures.At the same point, we realised that different accreted structures are likely occupying the same abundance space -as previously suggested by both Horta et al. (2021) and Rey et al. (2023).Only smallest measurement uncertainties avoid the blurring of sequences (Fig. 4).
• When tracing three sequences in the age-[Fe/H]-plane, we confirmed that their significant overlap appeared when these were plotted in terms of dynamical and chemical properties (Fig. 5).However, in the age-[Fe/H]-plane these sequences were clearly separated.We thus inspected the age-abundance relations in more detail (Sec.4).
• Plotting elemental abundances against age, proved fruitful, with very distinct sequences of accreted and in situ populations appearing in the simulation data (Figs.6 and 8) with up to 12-sigma separation (Tab.A1) at the youngest ages.
• Such sequences are not as clear in the observational plots of GALAH DR3 (Fig. 7), as a result of massive uncertainties that come with the probabilistic nature of estimating stellar ages through isochrone fitting.
• When testing the significance of separations of sequences in the age-[X/H] plane, we see the largest separations for the youngest accreted stars, even at larger observational noise (Fig. 9 and Tab.A2).We find that we loose the ability to tell sequences apart beyond 2-sigma for observational noise above 0.15 dex (Fig. 9) and age uncertainties above 20% (Fig. 10).
• With observational data below these uncertainties, we expect to see accreted and in-situ sequences separated (Fig. 11), as has been the case in the work by Xiang & Rix (2022).
• In this work, we have only analysed one Milky Way analogue with one clear in-situ sequence and one major as well as one minor accreted sequence.For the major accreted sequence, we find a stellar mass ratio of 1:3.0 at the time of last star formation in the accreted sequence (Fig. 14).
Future Research: These insights provide clear avenues for following up -both in the observational data and other simulations beyond the one Milky Way analogue of this study: • We can use the observed chemical evolution in both in-situ and accreted components -formed in systems with different masses and stellar mass distributions -to adjust the chemical yields and thus chemical evolution in simulations (see Fig. 2 for example for Ba as a neutron-capture element).
• Simulations are an excellent test bed to probe the efficiency of our observational tools to identify accreted stars.In particular, we plan to quantify the accuracy and completeness of previously used Gaussian Mixture Models (Das et al. 2020;Buder et al. 2022) with simulated [Al/Fe] vs. [Mg/Mn] plane in a follow-up paper.
• This study underscores the imperative for more precise and accurate age determinations.While better precision can be expected from additional information by asteroseismic missions (see e.g.Miglio et al. 2017;Mackereth et al. 2021), better accuracy can only be provided by more advanced stellar evolutionary models (see also Kim et al. 2002;Schuster et al. 2012) as an undeniable bottleneck of Galactic archaeology.
• We need to quantify the separation of sequences in other NI-HAO galaxies (see e.g.Lu et al. 2022;Buck et al. 2023) as well as other simulations, like VINTERGATAN (Renaud et al. 2021a,b;Agertz et al. 2021), HESTIA (Khoperskov et al. 2023b,c,d) and if possible beyond just age-[Fe/H] ( Khoperskov et al. 2023d).Constraining several of the simulation parameters, for example by up-or down-scaling merger significance to the point where they do or do not resemble our Milky Way properties is another exciting avenue along this line Rey et al. (2023).
• Following up the correlation and causality of the spread in the age-[Fe/H] relation among stars born after the major merger (Fig. 13) is an exciting avenue of further research.It will be important to compare the stellar, gas, baryonic, and total masses, inclination angles, merger timescales, strength and emergence of bars and other properties of accreted galaxies in simulations to the impacts of massive mergers on the cold gas metallicity gradients (Buck et al. 2023) as well as the correlation with observational plots that use age and/or chemistry.This will be complementary to the ongoing work in the kinematic/dynamic comparisons (e.g.Naidu et al. 2021;Khoperskov et al. 2023c) when reconstructing the complex pattern of accretion in the Milky Way (Naidu et al. 2020).
Fuelled by the revolutionary orbit information of the Gaia satellite (Gaia Collaboration et al. 2016Collaboration et al. , 2018Collaboration et al. , 2021Collaboration et al. , 2023) ) with its unprecedented astrometric (Lindegren et al. 2016(Lindegren et al. , 2018(Lindegren et al. , 2021) and spectroscopic information from both Gaia (Katz et al. 2019(Katz et al. , 2022) ) as well as ground-based surveys (e.g.Abdurro'uf et al. 2022;Buder et al. 2021;Zhao et al. 2012;Conroy et al. 2019), we are currently in the discovery phase of substructure in our Milky Way from observational data; with new structures being identified in diverse and creative ways (e.g.Naidu et al. 2021;Malhan & Rix 2024).This surge in discovery is undeniably thrilling, yet it is worth to keep in mind that our bigger goal remains to quantify the role of accretion on shaping our Galaxy.A major question continues to be, if accretion was responsibly for the change in star formation efficiency around the end of thick disk and onset of thin disk formation (Conroy et al. 2022) -or just happening at the same vivid time in our Galaxy's evolution (see the reviews by Helmi 2020;Deason & Belokurov 2024, for more extensive discussions).Our best way forward to answering this question is to connect the unique insights from resolved observations in the Milky Way with extragalactic observations and theoretical predictions from simulations (see e.     Table A1.Mean values and standard deviations of [X/H] in different age bins for the three age-abundance sequences of the Milky Way analogue as identified in Eqs.4-6 as well as the separation significance  12 between Sequence 1 and 2 as defined in Eq. 8. Values are calculated for the observable footprint without observational noise if more than 10 particles are within a bin.We do not show Ne and P, as they are commonly not measured by stellar surveys.Element Quantity 8.5-9.0 9.0-9.5 9.5-10.010.0-10.5 10.5-11.011.0-11.5 11.5-12.012.0-12.512.5-13.013.0-13.5 [

Figure 1 .
Figure 1.Density distribution of all particles (left panels) and those in a GALAH-like footprint (right panels) for particles of the NIHAO simulation g8.26e11 (color-coded by logarithmic density).Panels a) and d): iron vs. silicon abundances as tracers of contribution from SNIa and CCSN.Panels b) and e): abundance combination [Al/Fe] vs. [Mg/Mn] used to trace accreted stars in observations.Panels c) and f): age-[Fe/H] relations.The selection in panels a) and e) are drawn with a polygon with vertices at (-0.9,0.05),(-0.65,0.025),(-0.65,0.045),and (-0.9,0.07).A convex hull polygon of the selected stars is then drawn in other panels.

Figure 2 .
Figure 2. Logarithmic density distribution of elemental abundances [X/Fe] versus [Fe/H] for the ten elements overlapping between GALAH and NIHAO.First and third rows show observed data from GALAH while second and fourth rows show the simulated data for the same abundance ratios in a GALAH-like footprint.Brighter colors (towards yellow) indicate higher densities.For better visibility, plot ranges are adjusted to the data rather than to be equal among each panel.Dashed lines indicate Solar values.Simulated abundance trends for N, Ne, P, and S are shown in Fig. A1.

Figure 3 .
Figure 3.Comparison of the diagnostic plot of odd-Z elements [Al/Fe] or [Na/Fe] vs. the ratio of CCSN and SNIa elements [Mg/Mn] in simulation (top) and observation (bottom), with accreted stars being expected in the upper left of the distributions.Left panels a and c use the same axis limits, whereas panel b is a zoomed version of the simulation.Panel d complements the picture, as [Al/Fe] is hard to measure in GALAH (see missing data points in top left of panel c).

Figure 4 .
Figure 4. Logarithmic density distribution of particles in the [Al/Fe] vs. [Mg/Mn] plane when adding increasing levels of noise to the observational footprint of the NIHAO simulation.Gray dashed lines indicate the Solar value, whereas the black dashed line indicates the visual centre in [Al/Fe] of the accreted sequence at −0.07 dex.Panel titles indicate the added noise level, relating to a decrease of separation significance in [Al/Fe] down to Δ ∼ 1.8 for panel d).Color maps have been saturated at maximum values of 100 for all panels for increased contrast and better comparability.

Figure 5 .
Figure 5. Logarithmic density distributions of chronochemodynamic properties of the three separate sequences identified in the age-[Fe/H] relation (first column) -blue for Sequence 1 in the top row, red for Sequence 2 in the middle row, and purple for Sequence 3 in the bottom row.The second column shows the spatial density of stars in -.The third column shows the kinematic extend of stars in the Toomre diagram of  vs. √  2 +  2 .The fourth and fifth columns show how the chemical distribution in the [Fe/H] vs. [Si/Fe] and [Al/Fe] vs. [Mg/Mn] plane, respectively.

Figure 6 .
Figure 6.Age-Abundance trends of the simulated Milky Way analogue in the observation footprint.Columns show the different the elements C, O, Mg, Al, Mn, and Ba (from left column to right column), respectively.Rows show different notations of abundances, that is, logarithmic abundance [X/Fe] relative to iron (top row), logarithmic abundance [X/H] relative to hydrogen (middle row), and linear number density ratio  X/H (bottom row).

Figure 7 .
Figure 7. Age-Abundance trends of the Milky Way Solar neighbourhood as observed by GALAH DR3.. Columns show the different the elements C, O, Mg, Al, Mn, and Ba (from left column to right column), respectively.Rows show different notations of abundances, that is, logarithmic abundance [X/Fe] relative to iron (top row), logarithmic abundance [X/H] relative to hydrogen (middle row), and linear number density ratio  X/H (bottom row).Error bars in middle row indicate median uncertainties for stars below and above 7 Gyr.

Figure 8 .Figure 9 .
Figure 8. Histograms of [Mn/H] in 0.5 Gyr age bins for the three age-abundance sequences (Eqs.4-6).Colors are the same as in Fig. 5. Panel titles show age bins and total number of particles within each bin.Inset text with same colors show contribution (in percent) of each sequence to the age bins as well as 16/50/84th percentile.Values are only calculated if more than 100 particles of sequence are within a given age bin.

Figure 10 .
Figure 10.Explanation of our gedankenexperiment (panel a) to quantify the remaining separation significance  12 for larger age uncertainties for different elements.Panel a) shows the age-[Fe/H] plane as example from which we select the two sequences with increasingly larger age bins from 8.5 − 9 Gyr to 8.5 − 12.5 Gyr.Panels b) and c) show the resulting separation significance  12 for different elements when selecting sequences based on the increasing bins (recalculated as relative bin sizes that can be understood as relative age uncertainties).

Figure 11 .
Figure 11.Comparison of age-metallicity relation for different scatter realisations for NIHAO simulations.Panel a) shows NIHAO without scatter, panels b-e show realisations of the same plot with increasingly more scatter sampling representative of uncertainties quoted in pioneering surveys (see text for more discussion).Panel f shows the GALAH DR3 age-metallicity relationship for comparison.

Figure 12 .
Figure 12.Relative contributions of Sequence 2 (panel a) and Sequence 3 (panel b) across the galactic area  vs. .The observational footprint within the Milky Way is indicated with a black circle and the galactic plane with a grey dashed line.The larger coverage of  vs.  with bin sizes 1 and 0.5 kpc, respectively, is inspired by the extragalactic study by Martig et al. (2021, their Fig.16).

Figure 13 .
Figure 13.Tracing the chemistry of particles in different regions of the age-[Fe/H] relation for the whole Milky Way analogue.Panel a) shows the age-[Fe/H] relation as grey density plot with the five rectangles indicating the particles that are traced.Panel b) shows the [Fe/H] vs. [Mg/Fe] plane with the density contours (or scatter points in the case of the intermediate pink region) of the different traced particles.Panel c) shows the [Al/Fe] vs. [Mg/Mn] plane for the same particles.

Figure 14 .
Figure 14.Absolute (panel a) and relative (panel b) stellar mass contributions to the observed solar neighborhood as a function minimum stellar age.Panel a) show the stellar mass ratio of the major accreted sequence to the in-situ sequence.

Figure A1 .
Figure A1.Same as Fig. 2, but for the elements N, Ne, P, and S that were not measured by GALAH.

Figure A2 .
Figure A2.Same as Fig. 6, but for the elements N, Ne, P, S, V, Cr, and Co.

Figure A3 .
Figure A3.Same as Fig. 13 but for the whole Milky Way analogue.

Figure A4 .
Figure A4.Same as Fig.13but also tracing the youngest stars of the Milky Way analogue.Note that the traced groups are separated at the end of the merger (8.2 Gyr), rather than before it (9.0Gyr).

Table 1 . Median abundance [X/H] for particles of the Milky Way ana- logue with
4 < Age / Gyr < 5 in the observation footprint.These were used to shift the the abundance to agree with an expected Solar value.