The relationship between age, metallicity, and abundances for disk stars in a simulated Milky Way galaxy

Galaxy formation is a violent, chaotic process. Small galaxies are stripped of their stars and accreted (Searle & Zinn, 1978) by their larger counterparts, cold gas is accreted through filaments (e.g., Kereš et al. 2005), and structures such as spiral arms and bars form and disrupt stellar orbits (e.g., Roškar et al. 2008). Yet, in some ways, galaxy formation is also very orderly. Across populations, we observe that star formation proceeds from the central regions to the outskirts (e.g., Bezanson et al. 2009; Carrillo et al. 2020) as exhibited by abundance gradients with radius (Nidever et al., 2014; Hayden et al., 2015; Kaplan et al., 2016; Belfiore et al., 2017; Weinberg et al., 2018), and disk scale heights decrease as galaxies evolve (Wisnioski et al., 2015). The combination of these orderly and disorderly processes gave rise to our very own Milky Way.

With Galactic archaeology, and specifically through using stars as our rosetta stone, we can begin to disentangle the detailed formation history of our Galaxy. This is because stars exhibit atmospheric abundances that generally reflect the chemistry of the gas from which they formed, barring dredge-up processes that bring elements made in the core of stars onto the surface (Iben, 1965). Although a star’s orbital and kinematic properties can change as it migrates from where it was born, its chemical fingerprint or inventory is an intrinsic property that remains relatively unchanged throughout its lifetime. With every successive stellar birth and death, or mass loss from stellar winds, the gas from which newer stars will be formed is further enriched. Thereby, stars will have locked within them the chemical fingerprints inherited from their environment, as expressed in the time and place of their formation.

In the last decade, we have been able to peer into the detailed chemical abundances of different stellar populations in the Milky Way with large spectroscopic surveys. Large multi-object spectroscopic stellar surveys, such as the Apache Point Observatory Galactic Evolution Experiment (apogee, Majewski et al. 2017), the Galactic Archaeology with HERMES (galah, De Silva et al. 2015), the Gaia–eso survey (Gilmore et al., 2012), and the Large Sky Area Multi-Object Fibre Spectroscopic Telescope (LAMOST, Cui et al. 2012) have expanded our stellar census of the Galaxy. These surveys have enabled an extensive characterization of the Milky Way’s two-component disk, historically called thin and thick disks due to their different spatial properties, and later on understood in terms of the low-$\alpha$ and high-$\alpha$ disks, respectively (Bland-Hawthorn et al., 2019). Earlier work (e.g., Chiappini et al. 1997) have shown that these two populations have different formation timescales, and with these large Milky Way surveys, it has been further solidified that the two components are spatially, kinematically, and chemically distinct, with varying contributions as a function of their location in the Galaxy (Hayden et al., 2015). Many groups have posited different origins for this chemical bimodality including clumpy star formation (Clarke et al., 2019), radial migration (Sharma et al., 2020), or the presence of a (gas-rich) merger (Grand et al., 2020; Lian et al., 2020; Buck, 2020; Agertz et al., 2021). Following on from this example, chemical abundances, especially for a large number of stars, e.g., N $>\rm 10^{5}$, are key to understanding how the Milky Way built its disk. While chemistry informs which stars formed together, age ultimately ties these processes to the temporal evolution of the Galaxy.

There have been many advances in the realm of age derivation for individual stars within large surveys. This has largely been driven by the availability of a set of reference objects with precision ages. The Kepler mission, in particular (Borucki et al., 2010), has provided asteroseismic ages for red giant branch stars (RGB) (Pinsonneault et al., 2014, 2018). The Gaia satellite (Gaia Collaboration, 2018) has also measured precision parallaxes, allowing age derivation from isochrones around the main sequence turn off (MSTO). Leveraging the mass-dependent, dredge-up-induced [C/N] ratio at the stellar surface (Masseron & Gilmore, 2015; Martig et al., 2016), data-driven approaches that use high precision reference ages have enabled ages to be determined from large samples of spectra (e.g. Ness et al. 2016; Ho et al. 2017).

Ages have provided the means to understand nucleosynthesis over time, and link the chemical evolution to the assembly history of a galaxy. There is a known age-[Fe/H] relationship for the Milky Way disk, wherein stars decrease in [Fe/H] with older ages (e.g., Twarog 1980; Soubiran et al. 2008). However, whether there is a clear trend is debated in observational studies: previous works (Haywood et al., 2013; Bergemann et al., 2014; Rebassa-Mansergas et al., 2016) found that there is substantial scatter in the relation, with Nissen (2015) finding this lack of correlation existing over a large age interval (e.g., 8 Gyr). This same study, however, found that there is a tight correlation between the age and [X/Fe] of stars. It is important to further investigate these age-individual abundance, or age-[X/Fe], relations at a given [Fe/H] in the disk because the tightness of these relations illustrate that (1) by inverting the age-[X/Fe] relation, abundances can serve as chemical clocks providing ages for stars, and that (2) [Fe/H] and age capture majority of the crucial information about a star e.g., [X/Fe], orbits, and birth location (Bedell et al., 2018; Ness et al., 2019; Jofré et al., 2020; Hayden et al., 2020; Espinoza-Rojas et al., 2021; Casamiquela et al., 2021; Sharma et al., 2022; Ratcliffe et al., 2022). Additionally, comparisons between the slope in the [Fe/H]-[X/Fe] relation and the slope in the age-[X/Fe] relation appear to group elements accordingly into three distinct nucleosynthetic sites: core-collapse supernova (SNe II), white dwarf explosion (SNIa), and stellar winds (Sharma et al., 2022). Furthermore, the slope in the age-[X/Fe] relation is sensitive to the location across the [Mg/Fe] vs [Fe/H] plane (Yuxi et al., 2021; Horta et al., 2021), but the scatter around these relations is consistently small. Both the high resolution studies at [Fe/H] = 0 dex from Nissen (2015) and Bedell et al. (2018) as well as the larger samples exploring a wider range of metallicities (Ness et al., 2019; Sharma et al., 2022) have shown that the intrinsic scatter of stars around their age-[X/Fe] relations is on the order of $<$ 0.01-0.04 dex. The small intrinsic scatter in the age-[X/Fe] relations is also fairly insensitive to the current location of the stars in the disk (Ness et al., 2019; Yuxi et al., 2021), which is a natural consequence of radial redistribution or migration, whereby stars at any given radius, R, are from a wide range of initial radii (Frankel et al., 2018, 2019).

The growing literature on the age-[X/Fe] relations in observations could constrain the nucleosynthetic processes that produce the elements in the universe and the star formation history and interstellar medium conditions that gave rise to the observed age-[X/Fe] relations and scatters in the Milky Way. This is because stars not only remember their birth sites through their detailed chemistry, but for an established disk, these birth sites, at a given age, can be ordered in radius by sorting in metallicity. Likewise, for a given metallicity, sorting in age probes systematically different birth radii. The observed low intrinsic scatter around the age-[X/Fe] trends at fixed [Fe/H] indeed corroborate this picture, wherein stars are formed with element abundances governed by their birth radii and time of formation.

It is therefore timely that we understand what gives rise to this relationship in a self-consistent manner. That is, investigating this relationship where we know the galaxy formation history and we know the nucleosynthesic processes that together bring about the resulting chemistry. Here, we utilize cosmological simulations as a means to examine the age-[X/Fe] relations across the disk. Within these simulations, we have access to the ages and chemistry of stars with a direct knowledge of birth properties while unhindered by observational uncertainties. Critically, recent cosmological zoom-in simulations (e.g. Wetzel et al. 2016, Grand et al. 2016, Buck et al. 2019, Wetzel et al. 2022) now implement physically motivated processes (e.g. sub-grid turbulent metal diffusion, stellar feedback, chemical enrichment from SNe II, SNIa, and stellar winds) at high resolutions that can distinguish individual star forming regions.

Galactic archaeology studies with simulated Milky Way-like galaxies have proven themselves instructive in demonstrating physical processes that can give rise to observed trends. For example, cosmological simulations show that gas-rich mergers can give rise to the low-$\alpha$ sequence seen in the Milky Way disk (Buck, 2020; Agertz et al., 2021). Yu et al. (2021) used the Feedback in Realistic Environments 2 (FIRE-2) simulations (Hopkins et al., 2018; Wetzel et al., 2022), to demonstrate that the thick disk formation phase in these Milky Way-like galaxies coincides with the star formation changing from a bursty mode to a steady mode. Nikakhtar et al. (2021) similarly used FIRE-2 Milky Way-like galaxies and explored the different families of stars in the solar neighborhood, determined through a Gaussian mixture modeling of kinematics and [Fe/H]. They found simulated analogs to the real and observed stellar population distinctions in the Milky Way solar neighborhood, and subsequently related these to having different origins. Specifically, they found a thin disk component that is young, a halo component consistent with early and massive accretion events, and three thick disk components with heated orbits due to satellite interactions. Bellardini et al. (2021, 2022) also looked at [Fe/H] and abundance gradients in FIRE-2 Milky Way-like galaxies, where they found that the [Fe/H] and abundances change from being dominated by azimuthal variations to being dominated by radial variations at lookback times 7–7.5 Gyr ago. This age broadly coincides with an earlier study by Ma et al. (2017) as the transition from a chaotic, bursty mode of star formation to a calmer, stable disk in the FIRE-1 galaxy m12i. As illustrated by these studies, the star formation history, gas inflow and outflow, nucleosynthetic processes, and environments that these galaxies live in are known, and therefore aid in understanding the possible origins of observed stellar population characteristics and interpret chemo-dynamical signatures in the Milky Way.

Therefore, in this study, we similarly take advantage of a cosmological zoom-in simulation of a Milky Way-like galaxy to explore and establish the relationship of the age, metallicity, and individual element abundances of stars. We use the Ananke Gaia synthetic survey (Sanderson et al., 2020) to understand the observed tight age-[X/Fe] trend at a fixed [Fe/H]. Specifically, we aim to investigate this trend at different [Fe/H] and locations in the disk, and to show where the simulations reproduce this observed relationship, where it breaks down, and how this gap could be bridged. Ultimately, we aim to be able to put into context the physical processes in the simulated galaxy that give rise to the age-[X/Fe] trend it exhibits. In Section 2 we describe the simulation data we use for this work and the galaxy m12i. In Section 3 we discuss the observational data we use for comparison. In Section 4 we explore the abundance and location distributions of stars with similar ages. In Section 5 we show the age-[X/Fe] trends for the various elements that are tracked in the simulations at different [Fe/H] and galactocentric radii. In Section 6 we compare the intrinsic scatter in the age-[X/Fe] trends between observations and simulations. Lastly, in Section 7 we discuss the implications of our results and provide a summary of this work. We note that for the majority ¹¹1In Section 4, we analysed one aspect of birth location of stars, i.e., their distributions with age, which is not public data. of the analysis, we used the publicly available information in the Ananke Gaia synthetic survey.

We briefly introduce the three sets of observational data that we compare to regarding the [X/Fe] distributions, age-[X/Fe] relations, and the intrinsic scatter, $\rm\sigma_{int}$, around these trends in the simulations.

First, we use the galah survey data release 3 (DR3, Buder et al. 2021). This has an associated value added catalog from Sharma et al. (2022) that contains ages from the BSTEP code (Sharma et al., 2018) which makes use of stellar isochrones. In total, galah DR3 has abundances for 30 elements in 588,571 nearby stars in the Galactic disk. We derived a sample from this data set consisting of 102,785 stars having applied the following quality cuts: signal-to-noise ratio $>$40, flag_sp=0, flag_fe_h=0, and $4000<\mbox{$T_{\rm eff}$}<7000$ K. We further apply a cut in [$\alpha$/Fe] vs [Fe/H] shown in Figure 3 and as similarly done in previous observational studies (e.g., Ness et al. 2019; Bland-Hawthorn et al. 2019), to focus on just the low-$\alpha$ disk, leaving us with 81,230 stars. We have the elements C, O, Mg, Si, and Ca in common with the galah data and compare their [X/Fe] distributions for mono-age populations with those from the simulations (see Section 4). The study by Sharma et al. (2022) also provide comparison with age-[X/Fe] relations and the $\rm\sigma_{int}$ around them at varying [Fe/H].

Second, we compare to the age-[X/Fe] relations from the set of 79 solar twin stars from Bedell et al. (2018), where the ages were derived from the Yonsei-Yale isochrones (Demarque et al., 2004). This sample of solar twin stars generally have ages $<$8 Gyr. The authors report measurements for 30 element abundances, where they have C, O, Mg, Si, S, and Ca in common with this study. With regards to the [X/Fe]-age trends, the fits were derived through minimizing the orthogonal distance between the data and the model, weighted by the observational uncertainties (Hogg et al., 2010).

Lastly, we use the measurements from the set of $\approx$ 15,000 red clump field stars from Ness et al. (2019) in the low-$\alpha$ disk, where the authors utilize the apogee survey’s 14th data release (Majewski et al., 2017). This work looked into the age-[X/Fe] relations in 19 different elements, wherein C, N, O, Mg, Si, S, and Ca are similarly measured as in the simulations. The ages were derived through data-driven modeling with The Cannon (Ness et al., 2015) using the APOKASC catalogue (Pinsonneault et al., 2018) as a training set which contains ages derived from asteroseismology. This sample of giant stars have ages $<$10 Gyr. In deriving the intrinsic scatter in the [X/Fe]-age trends in their study, the effect of the age uncertainties were taken into account via Monte Carlo sampling. The error in age can shift stars horizontally in the age-[X/Fe] trend. For flat trends, accounting for the age error has minimal effect but for steep trends, this essentially flattens the relationship. This comparison observational work found that the effect of the age uncertainties changes the intrinsic scatter by $<$0.01 dex and was therefore negligible. Likewise, the abundance error was accounted for in deriving the $\rm\sigma_{int}$ as the measured scatter from the best fit line will be a combination of the $\rm\sigma_{int}$ and the error in the abundance measurement, $\rm\sigma_{abund}^{2}$. That is: $\rm\sigma_{total}^{2}=\sigma_{abund}^{2}+\sigma_{int}^{2}$.

A starting point to understanding the present day distributions of abundances across the disk (see Section 5) is to look at the birth radius distributions for different age population. In Figure 4, we show the kernel density estimation (KDE) for the birth radius, $\rm R_{birth}$, of stars centered on different age bins from ages = 1 to 9 Gyr. Because the number of stars within each age bin is different, we normalize the KDEs such that the distributions are on the same scale for easier comparison. We used the smallest possible bin in age such that the measured dispersion across the [Fe/H] vs $\rm R_{birth}$ trend is not induced by the age binning, giving age bin sizes of (0.063, 0.125, 0.063, 0.125, 0.250) Gyr for the (1, 3, 5, 7, 9) Gyr samples, respectively. For more details on this process, we refer to Appendix B. We applied this same age binning in comparing the KDEs at different ages for [Fe/H] and [X/Fe], where X=C, N, O, Mg, Si, S, Ca shown in Figure 5.

The $\rm R_{birth}$ distributions peak from the outskirts of the galaxy to the central regions, from the youngest to the oldest stars, respectively. However, this is not observed for the 9 Gyr stellar population, which is more spread out across the disk. This is because at higher redshifts (i.e. lookback time $>$ 8 Gyr), the star formation in m12i was violent and bursty until a final minor merger finishes at lookback times of $\sim$6 Gyr, after which a more stable and calm disk appears and persists (Ma et al., 2017).

Next, we focus on the distributions of element abundances as a function of age as shown in Figure 5.

Metallicity: The [Fe/H] KDEs follow an increasing trend with decreasing age, though the 1 Gyr and 3 Gyr age bins have similar [Fe/H] distributions, and the width of the distribution also varies with different age bins. The similarity in the [Fe/H] distributions for the younger stars is in line with the little evolution in [Fe/H] and the similar radial abundance gradients at these ages found by Bellardini et al. (2022) in their sample of Milky Way-like galaxies in FIRE-2. Qualitatively, the [Fe/H] distribution widths seem to be related to the $\rm R_{birth}$ distribution widths shown in Figure 4, although we caution that this could be due to our age-binning, which was derived from the [Fe/H] vs $\rm R_{birth}$ trend at different ages. For the 9 Gyr bin, there is star formation everywhere in the galaxy, which induces a wide distribution in [Fe/H]. There is also a large tail towards lower metallicities reflecting the chemical enrichment at earlier times. Additionally, the large tail is affected by the merger with another smaller system, seen as a slight bump in the $\rm R_{birth}$ distribution at 20 kpc in Figure 4. For the 7 Gyr and 5 Gyr age bins, their [Fe/H] distributions are narrower, following similarly narrow and more-peaked distributions in $\rm R_{birth}$. Lastly, for the 3 Gyr and 1 Gyr stellar populations, the [Fe/H] distributions are similar, peak above solar values, and have broader distributions, reflecting the star formation that is happening in a wider range of $\rm R_{birth}$ across the disk. In observations, there is a known age-metallicity relationship in Milky Way disk stars (Soubiran et al., 2008), that shows decreasing [Fe/H] with older ages. On the other hand, many studies (Bergemann et al., 2014; Rebassa-Mansergas et al., 2016) have also found that there is no clear trend in this relationship. In a sense, the [Fe/H] distributions in this work agree with both scenarios, with the important distinction being the age at which this holds true. In this work, we find that the 1 and 3 Gyr age bins show little change in their respective [Fe/H] distributions, while increasing in age from the 5 Gyr bin shows progressively more metal-poor distributions. Bellardini et al. (2022) shows this more clearly for their sample of Milky Way-like galaxies in FIRE-2 in looking at the age-metallicity relationship in the simulations (see their Figure 1). The [Fe/H] shows no trend with age for the youngest stars up to $\sim$4 Gyr, after which there is a slight decrease in [Fe/H] with increasing age up to 7 Gyr, and a more drastic decrease in [Fe/H] at much older ages.

Looking at the trends in [X/Fe] distributions with age in Figure 5, we find some elements follow expectations from observations while others diverge. We discuss these [X/Fe] distribution trends in order of atomic mass, from the light elements to the $\alpha$-elements. We also list the summary statistics of abundance distributions in Table 1 and where available, we included summary statistics from observations of the Milky Way low-$\alpha$ disk taken from galah DR3 written in brackets. We note however that the age bins used for the observations are much bigger (e.g., $\pm$0.8 Gyr) to take into account age uncertainties. Due to a relatively smaller sample, we also did not include a 1 Gyr bin from the observations.

Carbon: The carbon abundance, [C/Fe], increases in peak [C/Fe] value with younger populations in the simulation. In this work, we have shown that stars in different age bins have roughly different [Fe/H] distributions in the -1 $<$ [Fe/H] $<$ 0 dex range; with this relation, we find that the [C/Fe] distributions in the simulations are contrary to what is found in observations, wherein [C/Fe] decreases with increasing metallicity at [Fe/H]$>$-1 dex (e.g. Bensby & Feltzing 2006). This is also supported by our comparison to the galah data as listed in Table 1, where the mean [C/Fe] increases with older ages in the observations but decreases in the simulations. C is mostly produced in massive stars, followed by low-mass AGB stars (Kobayashi et al., 2020) in the Galaxy. In the simulations, there is a larger contribution of C coming from stellar winds, which is one possible source of discrepancy. In observations, the atmospheric abundance of C is also affected by dredge-up in convective stars. Throughout the life of a star, its C+N+O abundances remain relatively constant. However, the slowest reaction in the CNO cycle produces a net increase in $\rm{}^{14}N$, decrease in $\rm{}^{12}C$, and slight decrease in $\rm{}^{16}O$ (Masseron & Gilmore, 2015; Martig et al., 2016). This is especially seen with RGB stars that have deep convective zones which take materials formed in the core onto the surface of the star, changing the atmospheric abundance. However, in the simulations, we have yields (not necessarily atmospheric abundance) which are IMF-averaged and all stars generated from a star particle inherit the same chemistry. Therefore, although this process is prescribed on the population level in the simulations (i.e., from van den Hoek & Groenewegen 1997, Marigo 2001, Izzard et al. 2004), the abundances for individual stars are more nuanced than can be achieved by this procedure.

Nitrogen: The nitrogen abundance, [N/Fe] distribution distinctly changes with the stellar population age, from broader with lower [N/Fe] abundances at older ages and narrower with higher [N/Fe] abundances at younger ages. The [N/Fe] distributions at different ages also overlap the least among the elements considered in this work, which is a result of its yield’s metallicity-dependence from SNe II as implemented in the simulations (Hopkins et al., 2018). In FIRE-2, only the elements N and O have yields from progenitors with metallicity-dependence, which is linear with progenitor [Fe/H] up to $Z/Z_{\odot}$ = 1.65. N is also intricately linked to C through the CNO cycle. However, compared to C, N is mainly produced in intermediate-mass AGB stars (Kobayashi et al., 2011). Like C, the effects of dredge-up on the N abundance is applied in the simulation yields, although IMF-averaged. Works using large spectroscopic survey data like apogee (e.g. Hayes et al. 2018) find a decrease in [N/Fe] with increasing [Fe/H] up to [Fe/H]$\sim$-0.5 dex, after which [N/Fe] increases with increasing [Fe/H], similar to what we find in the simulations.

Oxygen: The oxygen abundance, [O/Fe], in Ananke shows the smallest evolution across time among the elements in this study i.e., the [O/Fe] distributions show the most insignificant differences for stars of different ages. The [O/Fe] distributions are wide at older ages and become narrow at younger ages, but, with each distribution overlapping substantially. In the simulations, most of the O is produced via SNe II which is expected for this hydrostatic $\alpha$-element. However, in FIRE-2, there is also non-negligible contribution from stellar winds with a dependence on the progenitor metallicity. The stellar winds contribute a comparable amount of O to the ISM as SNe II at 3 $\times$ solar [Fe/H]. This trend in the simulations is contrary to what is found in observations, where [O/Fe] is seen to decrease with increasing [Fe/H] (e.g. Bensby et al. 2004), and loosely, age, as is expected of an $\alpha$-element. This increase in abundance with older age is indeed what we also find for the [O/Fe] distribution trends in galah. The significant O production from stellar winds at higher metallicities is a possible source of discrepancy.

$\alpha$-elements (Magnesium, Silicon, Sulfur, Calcium): Lastly, we consider the other $\alpha$-element abundance distributions i.e. Mg, Si, S, and Ca. Although these elements have distributions with different mean values, the characteristics of the distributions for the different age bins are similar: broad at older ages with higher abundances, evolving to narrow peak at younger ages with lower abundances. This behavior we see for the mean [X/Fe] is expected for $\alpha$-elements, which are mainly produced through SNe II with shorter timescales, and diluted via SNIa that have delayed and longer timescales, therefore decreasing the $\alpha$-abundance with younger generations of stars. This is also what we observe in the galah data, aside from Ca, which is intriguingly constant with age in observations (see also discussion in Section 5), but increasing with older ages in the simulations like the other $\alpha$-elements. Mg, like O, is a hydrostatic $\alpha$-element produced during the shell burning phase in massive stars. Si, S, and Ca are explosive $\alpha$-elements that also have non-negligible contributions from SNIa. For the explosive $\alpha$-elements, the peaks in the [X/Fe] distributions between the oldest and youngest stellar populations are separated by only $<$0.1 dex in the simulations. On the other hand, the same mono-age populations are more widely separated for Mg that has negligible contributions from SNIa. In the simulations, Mg and Si are produced in similar amounts from SNe II, with S yields being 1/3rd and Ca yields being 1/25th of the Mg yields. Among the $\alpha$-elements, SNIa contributions are highest for Si and S, and lower by a factor of 13 and 20 for Ca and Mg, respectively. We also note a tail in the distributions of these elements towards lower [X/Fe], at all age bins, and most pronounced at 1 Gyr.

In this section, we found that mono-age populations in the thin disk of m12i within $\rm R<25kpc$ trace birth locations (Figure 4) that overlap quite broadly but peak at different $\rm R_{birth}$ and exhibit [Fe/H] and [X/Fe] distributions that also show distinct trends (Figure 5). Interestingly, the [Fe/H] trends mirror $\rm R_{birth}$ trends, where a wider $\rm R_{birth}$ distribution also corresponds to a wider [Fe/H] distribution; together, age and [Fe/H] trace $\rm R_{birth}$ in the disk of m12i.

Irrespective of the [Fe/H] distributions, however, the width of the [X/Fe] distributions increases at older stellar ages. The oldest stars (ages $>$9 Gyr) were born from less homogeneous gas, when star formation was happening more chaotically in the galaxy as seen in previous works (e.g. Ma et al. 2017, Agertz et al. 2021, Bellardini et al. 2021). This is further supported by the larger azimuthal abundance scatter at higher lookback times for similar Milky Way-like galaxies in FIRE-2 (Bellardini et al., 2022). This subsequently produces wider abundance distributions. The increasing kurtosis in [X/Fe] with younger stellar ages, best seen in the elements N and O, indicates that the chemical enrichment in the disk becomes more homogeneous over time. At the same time, the skew in the [X/Fe] distributions becomes more pronounced at younger ages, which could be explained by where the stars were born. Focusing on the 1 Gyr population, part of the sample was born at smaller radii, from gas that has had further dilution from SNIa (i.e., greater overlap with older populations as seen in Figure 4), therefore skewing the [X/Fe] distribution more negatively. However, these populations can be disentangled with [Fe/H], which would tell us which stars were born more centrally (i.e. higher [Fe/H]) or in the outskirts (i.e. lower [Fe/H]). Therefore, age is indeed intricately linked to the [Fe/H], [X/Fe], and $\rm R_{birth}$ of a star.

In this section, we also showed that the trends in the [X/Fe] distributions with age for the elements N, Mg, Si, and S in the simulations agree with expectations from observations, while those for the elements of C, O, and Ca do not. These similarities and differences have repercussions in the resulting age-[X/Fe] trends that we discuss in the following section.

We have seen in the prior section the spatial trends of stellar element abundances, that are imprinted at birth. However, stars of different ages will have since moved from their birth location, so any initial trends will not be trivially maintained. As shown in Figure 6, there is a smoother distribution of current stellar location, $\rm R_{now}$ compared to $\rm R_{birth}$, as blurring and churning (i.e. radial migration) processes wash out the original clumpiness that stars were born into (e.g., Frankel et al. 2020). However, both $\rm R_{now}$ and $\rm R_{birth}$ still show distinct trends at different age groups with younger stars in the outskirts and older stars more centrally-located.

In this Section, we examine what the current distribution of element abundances at different locations and [Fe/H] as a function of age looks like, as a result of the combination of birth properties and migration. We note that from within the cuts outlined in Section 2.2, we pivot to two fixed [Fe/H] values, of [Fe/H] = 0 dex and [Fe/H] = $-0.25$ dex, for the purpose of demonstrating the age–[X/Fe] trends for a high and a low [Fe/H] sample, respectively.

We take a deeper look at the $\rm\sigma_{int}$ in the simulation around the age-[X/Fe] relations, compared to Milky Way, as well as at different [Fe/H] and locations. We aim to see if this measurement of the scatter around the age-[X/Fe] relations can be used as a diagnostic and/or as a metric to capture the chemical enrichment in the simulated Milky Way galaxy m12i, and to relate this quantity to what is happening in the galaxy at that point in time.

The stars in the Milky Way have been shown to exhibit tight age-[X/Fe] relations at a given [Fe/H], such that by knowing a star’s age and [Fe/H], one can predict its abundance to a precision of $\sim$0.02 dex. Indeed, this has been studied in greater detail and in many axes e.g., high vs low-$\alpha$ disk, high vs low [Fe/H], MSTO vs RGB stars, thanks to large spectroscopic survey data (Ness et al., 2019; Yuxi et al., 2021; Sharma et al., 2022).

The main goal of this work is to explore this observed age-[X/Fe] relations in cosmological zoom-in simulations that have been shown to reproduce Milky Way results from Gaia DR2 (Sanderson et al., 2020). One advantage to doing this is that we have absolute knowledge of when and where stars formed in the simulations. In contrast, ages in observations are harder to derive, the approach for deriving them is dependent on the type of star, and birth locations are unknown. Our main focus has been to determine if the current state-of-the-art simulations that trace the chemical enrichment of galaxies, can successfully reproduce the observed age-[X/Fe] trends in the Milky Way, and where they fall short and disagree. However, we can go beyond this and even explore what is physically happening in the simulated galaxy, and investigate where stars of a given age or chemistry are born.

Therefore, in understanding where the age-[X/Fe] trends in the simulations come from, we have also come to establish the intricate relationship of age, metallicity, individual element abundances, and location of stars in this simulated Milky Way galaxy, m12i. We have shown that aside from the chemical enrichment in the galaxy being a product of the prescribed nucleosynthesis (see Section 4), the individual element abundances, and specifically their scatter around the age-[X/Fe] relations, $\rm\sigma_{int}$, generally reflect m12i’s star formation history. Focusing on the solar [Fe/H] sample, at all distances, we find that the running $\rm\sigma_{int}$ across age is constant at young ages, up to about $\sim$4 Gyr, and then increases to higher $\rm\sigma_{int}$ values at older ages. The detailed study of m12i by Ma et al. (2017) find that the disk changes from a bursty and clumpy star formation mode to a stable star formation mode at look back times of 6 Gyr. This galaxy was also included in the more recent analysis by Yu et al. (2021) where they found that m12i changes to a near constant state of star formation at $\sim$3 Gyr. Indeed, for stars that are $<$3 Gyr, the $\rm\sigma_{int}$ is constant, for 3$<$age$<$6 Gyr, the $\rm\sigma_{int}$ increases, and at $>$6 Gyr, the $\rm\sigma_{int}$ is the largest. We would like to point out however, that we mainly see this behavior for the solar [Fe/H] sample.

To make sense of the observed age-[X/Fe] trends, we have explored their distributions at different ages in Section 4. Importantly, we also explored where these stars are born and where they are currently as shown in Figures 4 and 6. Here we find that stars in the oldest age bin (9 Gyr) are spread all across the disk within R $<$ 10 kpc, extending farther than stars that are younger in the 5 Gyr and 7 Gyr age bins. This might seem counter-intuitive in the inside-out formation scenario, where star formation starts in the central regions and proceeds outwards. However, for these oldest stars, the galaxy that eventually forms into m12i is composed of multiple smaller galaxies (see Figure 1 in Ma et al. 2017), which gives rise to the more extended and spread out $\rm R_{birth}$ distribution that we see. Indeed, Santistevan et al. (2020) found that the Milky Way/M31-mass galaxies in the FIRE simulations are made up of a hundred distinct dwarf galaxies before z$\sim$2. Therefore, the stars at the oldest age bin show the largest spread in the abundance distributions (as similarly seen in the analysis from Bellardini et al. 2022 of the stellar abundances of Milky Way/M31-mass galaxies), which reflects the chemistry of this clumpier and more chaotic star formation, and with the [Fe/H] KDE showing a long tail towards lower values. Curiously, the stars in the 3 Gyr bin have a similar $\rm R_{birth}$ (as well as current R) distribution to the 9 Gyr age bin, but its abundance distributions are narrower. This could be explained however by the existence of a rotationally-supported and stable disk that at 3 Gyr, has stars with abundances that are more azimuthally homogeneous at a given radius (Bellardini et al., 2022).

We also note that stars with different ages have significant overlap in their distributions of [X/Fe], $\rm R_{birth}$, and current R, but that the peak of those distributions are distinct from each other. For younger stars, the peak of the radial distributions of stars is at larger R, and the abundance distributions exhibit progressive enrichment. There is an overlap in the [X/Fe] KDEs, which is manifested in the gentle slopes in the age-[X/Fe] trends. The element N, whose yields from SNe II have a progenitor metallicity-dependence, has the most distinct [N/Fe] distributions as a function of age (i.e., [N/Fe] distributions overlap the least at varying age bins), and therefore also the steepest age-[X/Fe] trend. On the other hand, O, whose yields from stellar winds have a progenitor metallicity-dependence, has abundance KDEs that overlap for different ages, and therefore exhibits the flattest trend with age. We therefore find that for the simulated Milky Way-like galaxy m12i, the direction of the age-[X/Fe] trends is a reflection of the chemical evolution prescriptions in the simulations, and the scatter around these relations is a reflection of the mode of star formation, with higher scatter at earlier times when star formation was more bursty and chaotic, and lower scatter at later times when the disk is more stable. Meanwhile, the magnitude of the slopes in the age-[X/Fe] relations is a convolution of the effects from the chemical evolution prescription as well as the star formation history. With the next generation of spectroscopic surveys dedicated to densely mapping stars in the Milky Way such as SDSS-V (Kollmeier et al., 2017), WEAVE (Bonifacio et al., 2016), and 4MOST (de Jong et al., 2019), as well as novel ways to derive ages for these stars, we can use the $\rm\sigma_{int}$ in the age-[X/Fe] trend as another diagnostic in understanding the evolution of the Milky Way.

In Section 6 we compared the slopes of and $\rm\sigma_{int}$ around the age-[X/Fe] relations from this work to those derived in observations. Specifically, in Figure 11 we show the $\rm\sigma_{int}$ in the solar neighborhood at solar metallicity for the stars in Ananke, solar twins (Bedell et al., 2018), and RGB stars (Ness et al., 2019). We find that the $\rm\sigma_{int}$ for m12i and the Milky Way are similar, and in fact, some elements such as Mg and Si have near-identical $\rm\sigma_{int}$ across all three samples. The elements C, O, and Ca, however, show the opposite trends in their slopes with age in the simulations compared to observations. The age-[C/Fe] relation increases with older ages in observations and has been noted to even have a steeper slope with age than $\alpha$ elements Mg and Si (Nissen, 2015). In contrast, we find that in m12i, the age-[C/Fe] relations have a negative slope; decreasing with increasing age. We expect that this is because of the substantial contribution from stellar winds at later times prescribed in the simulations. The age-[O/Fe] relations show a positive slope in observations, as is expected for an $\alpha$-element. However, we find that the O abundance is flat, or decreasing with older ages in the simulations. This is due to the large O yields from stellar winds at higher metallicities as prescribed in the simulations. Lastly, observations find that the [Ca/Fe] trend with age is flat or decreasing with older ages, contrary to the trend for the other $\alpha$ elements (Nissen, 2015; Bedell et al., 2018; Ness et al., 2019; Sharma et al., 2022) and to what we find in this study, though this discrepancy could potentially be alleviated by the inclusion of other astrophysical sources (e.g., Mulchaey et al. 2014). The chemical enrichment in this cosmological zoom-in simulation of a Milky Way-like galaxy and observational comparisons highlight that our knowledge of how elements form and evolve through time is incomplete. Continued efforts for Galactic chemical evolution modeling (e.g. Kobayashi et al. 2020) in tandem with a more flexible way of propagating chemical enrichment are therefore important for these simulated galaxies to be an even more realistic laboratory for Galactic archaeology. Nonetheless, it is truly noteworthy that m12i, which was evolved in a large cosmological volume (later on zoomed-in which the Ananke synthetic survey was made from), and was matched as a Milky Way galaxy based only on its mass and isolation, reproduce comparable age-[X/Fe] trends to the Milky Way.

Moving forward, there are many avenues to delve deeper into Galactic archaeology in simulations, especially in the Ananke catalog, to link observed Milky Way properties to a plausible galaxy formation history. For example, it is worth looking more into the $\rm\sigma_{int}$ for the different Milky Way galaxies in Ananke that have masses 1-2 $\times 10^{12}\rm M_{\odot}$ but have varying morphologies and star formation histories. This would be important in contextualizing what brings about the observed $\rm\sigma_{int}$ in the Milky Way by comparing to the derived $\rm\sigma_{int}$ in simulated Milky Way galaxies with different properties. In line with this, we have shown in Section 6 how the $\rm\sigma_{int}$ for observations and simulations compare in the solar neighborhood, with the former calculated from galah abundances of Milky Way stars and the latter based on this analysis. With future surveys that will observe stars in the Milky Way more extensively (e.g. Kollmeier et al. 2017), we will be able to map the $\rm\sigma_{int}$ in the Galaxy in its different regions, that have, in turn, their distinct stellar populations to aid in a more holistic picture of its formation history. Bellardini et al. (2021) used simulated Milky Way galaxies to show that the abundance scatter in gas is dominated by radial trends at lookback times $<$6.9 Gyr, and follow-up work show that this transition lookback time is $<$7.5 Gyr for stars (Bellardini et al., 2022), beyond which the abundance scatter is dominated by azimuthal inhomogeneity. It is therefore worthwhile to also explore the age-[X/Fe] trends from different solar viewpoints in the same galaxy to ultimately understand if the age-[X/Fe] relationship at a given metallicity is dominated by radial or azimuthal trends.
The future of the field is bright with large spectroscopic surveys in tandem with novel ways of deriving ages, as well as updates on hydrodynamical simulations to produce more realistic star formation in galaxies. This will allow the community to tackle how our Galaxy formed more expansively and from different perspectives, bridging our knowledge and the gap between observations and simulations, and with stars, specifically their chemical fingerprint, as the main tool. In this work, we give but a glimpse of this future, outlining how the abundance $\rm\sigma_{int}$ can be a diagnostic for unraveling the evolution of the Milky Way. In the end, we are still left with many questions to investigate (e.g. galactic chemical evolution modeling, azimuthal vs radial abundance trends, flexible chemical prescriptions) both from the observational and simulation point of view, and we look forward to answering them in the future, in this era of big data in Galactic archaeology.

AC is grateful to the Flatiron Institute Center for Computational Astrophysics Pre-Doctoral Program through which this project was made possible. AC also acknowledges support from the Science and Technology Facilities Council (STFC) [grant number ST/T000244/1] and thanks the Large Synoptic Survey Telescope Corporation (LSSTC) Data Science Fellowship Program, which is funded by LSSTC, NSF Cybertraining Grant 1829740, the Brinson Foundation, and the Moore Foundation. MKN is in part supported by a Sloan Research Fellowship. AW received support from: NSF via CAREER award AST-2045928 and grant AST-2107772; NASA ATP grant 80NSSC20K0513; HST grants AR-15809, GO-15902, GO-16273 from STScI. KH and AC acknowledge support from the National Science Foundation grant AST-1907417 and AST-2108736 and from the Wootton Center for Astrophysical Plasma Properties funded under the United States Department of Energy collaborative agreement DE-NA0003843. This work was performed in part at Aspen Center for Physics, which is supported by National Science Foundation grant PHY-1607611. This work was performed in part at the Simons Foundation Flatiron Institute’s Center for Computational Astrophysics during KH’s tenure as an IDEA Fellow. RES acknowledges support from the Research Corporation through the Scialog Fellows program on Time Domain Astronomy, from NSF grant AST-2007232, from NASA grant 19-ATP19-0068, and from HST-AR-15809 from the Space Telescope Science Institute (STScI), which is operated by AURA, Inc., under NASA contract NAS5-26555.

Simulations in this project and the mock catalogs generated from them were run on the Caltech compute cluster “Wheeler,” allocations from XSEDE TG-AST130039 and PRAC NSF.1713353 supported by the NSF, NASA HEC SMD-16-7592, and the High Performance Computing at Los Alamos National Lab., and analyzed using computing resources supported by the Scientific Computing Core at the Flatiron Institute. Further simulations were run using Early Science Allocation 1923870 on Frontera, which is made possible by National Science Foundation award OAC-1818253. This work used additional computational resources of the University of Texas at Austin and TACC, the NASA Advanced Supercomputing (NAS) Division and the NASA Center for Climate Simulation (NCCS), and the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation grant number OCI-1053575. Last but not the least, huge thanks to Jason Hunt who generously hosted the data as we were closing in on the finish line.