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

Observations of the Milky Way's low-$\alpha$ disk show that at fixed metallicity, [Fe/H], several element abundance, [X/Fe], correlate with age, with unique slopes and small scatters around the age-[X/Fe] relations. In this study, we turn to simulations to explore the age-[X/Fe] relations for the elements C, N, O, Mg, Si, S, and Ca that are traced in a FIRE-2 cosmological zoom-in simulation of a Milky Way-like galaxy, m12i, and understand what physical conditions give rise to the observed age-[X/Fe] trends. We first explore the distributions of mono-age populations in their birth and current locations, [Fe/H], and [X/Fe], and find evidence for inside-out radial growth for stars with ages<7 Gyr. We then examine the age-[X/Fe] relations across m12i's disk and find that the direction of the trends agree with observations, apart from C, O, and Ca, with remarkably small intrinsic scatters, $\sigma_{int}$ (0.01-0.04 dex). This $\sigma_{int}$ measured in the simulations is also metallicity-dependent, with $\sigma_{int}$ $\approx$ 0.025 dex at [Fe/H]=-0.25 dex versus $\sigma_{int}$ $\approx$ 0.015 dex at [Fe/H]=0 dex, and a similar metallicity dependence is seen in the GALAH survey for the elements in common. Additionally, we find that $\sigma_{int}$ is higher in the inner galaxy, where stars are older and formed in less chemically-homogeneous environments. The age-[X/Fe] relations and the small scatter around them indicate that simulations capture similar chemical enrichment variance as observed in the Milky Way, arising from stars sharing similar element abundances at a given birth place and time.


INTRODUCTION
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 andreia.carrillo@durham.ac.uk 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-α and high-α 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 , 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 > 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(Pinsonneault et al. , 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 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 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(Frankel et al. , 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) 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-α 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. (2021Bellardini et al. ( , 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 char-acteristics 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 Waylike 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 1 of the analysis, we used the publicly available information in the Ananke Gaia synthetic survey.

Ananke: Gaia Synthetic Surveys
We take advantage of the Gaia synthetic surveys produced with the Ananke framework (Sanderson et al. 2020) which generates synthetic phase-space surveys from baryonic simulations. These synthetic surveys are based on the Latte suite of zoom-in simulations of Milky Way-like galaxies (Wetzel et al. 2016) utilizing FIRE-2 physics (Hopkins et al. 2018), that features state-of-theart implementations of hydrodynamics, radiative cooling and heating, star formation, and stellar feedback.
In this simulation suite, stars are formed from selfgravitating molecular gas that is self-shielding (following Krumholz & Gnedin 2011), dense (n > 1000 cm −3 ), cold (T< 10,000 K), and Jeans-unstable-prescriptions that produce clustered stellar populations naturally. A star particle is essentially a single stellar population with an initial mass of 7070 M but with mass loss from stellar winds, it reduces to a typical mass of ∼5000 M at z = 0.
The stellar evolution of star particles is produced from STARBURST99 (Leitherer et al. 1999) with a Kroupa (2001) initial mass function (IMF). Chemical enrichment occurs through three main nucleosynthetic sources: SNe II, SNIa, and OB/asymptotic giant branch (AGB) stellar winds. The yields for SNe II are from Nomoto et al. (2006), for SNIa from Iwamoto et al. (1999), and OB/AGB winds from Wiersma et al. (2009) compiling values from van den Hoek & Groenewegen (1997), Marigo (2001), and Izzard et al. (2004). SNIa rates are taken from Mannucci et al. (2006) including both prompt and delayed populations. We note that the yields in the simulations are IMF-averaged (Hopkins et al. 2018). These stellar and chemical evolution prescriptions (listed in Appendix A) affect the present-day "observed" abundances in the synthetic survey which is central to our analysis.
Sub-resolution gas metal diffusion is included, which leads to abundance distributions that better match observations (Su et al. 2017;Hopkins et al. 2018;Escala et al. 2018). For a more detailed exploration of how the diffusion coefficient affects abundance scatters in Milky Way-like galaxies, see Bellardini et al. (2021), Appendix C.
The simulations track the total mass of each element for each star particle, but observationally, the abundance measurement is in terms of [X/Fe], where [X/Fe] = [X/H] -[Fe/H] and the bracket,"[ ]", notation denotes relative abundance with respect to Solar abundance values. To put the abundances from the simulations in the same scale, Asplund et al. (2009) Solar values were used to calculate where X is the element, m X is the mass of the element in the star particle, and m X, is the Solar value for that element. Three galaxies were taken from Latte, labeled m12f, m12i, and m12m, for the generation of the Gaia synthetic surveys viewed from three different "Solar" view points per galaxy (resulting in nine synthetic surveys in total), where a "Solar" view point (Local Standard of Rest, LSR) is defined as 8.2 kpc from the galactic center. These Milky Way-like galaxies were selected based only on their total mass (M 200m 2 = 1-2 ×10 12 M ) and environment (isolated, i.e., has no nearby dark matter halo with similar mass up to 5 ×R 200 3 ) and therefore have a variety of morphologies, satellite populations, and star formation histories. Each catalog has self-consistent dust extinction implemented, based on scaling of the gas density of the simulated galaxies. For a proof of concept in exploring age-[X/Fe] trends in simulations, and for simplicity, we focus our analysis on one viewpoint (LSR 2) of the galaxy m12i which is also publicly available 4 . The Galactic map of m12i is shown in Figure  1.
For the synthetic surveys, the single stellar populations are broken into individual synthetic stars that are sampled similarly from a Kroupa (2001) IMF and mapped to a grid of isochrones in age and [Fe/H] over a set of model isochrones (PARSEC; Bressan et al. 2012,Marigo et al. 2013,Marigo et al. 2017. We note that the isochrones used for generating the synthetic surveys are different from the isochrones used for the stellar evolution in the simulations. Specifically, they differ in their treatment of high-mass star evolution where the simulations use the Geneva tracks (Lejeune et al. 1997) that have enhanced mass loss from rotation while PARSEC includes thermally-pulsating AGB stars. However, Sanderson et al. (2020) prescribed massindependent mass-loss rate in sampling the number of stars from a star particle, alleviating some of the differences introduced from the mismatch in isochrones.
These individual stars then inherit the age, [Fe/H], and [X/H] of the parent star particle, while their 6D phase-space information is sampled over a onedimensional kernel in each position and velocity axis centered on the parent star particle. We use the true values for the "observed" properties (e.g. distance, age, and chemistry) of stars in the synthetic survey because the essence of our analysis is to understand the intrinsic scatter in the age-[X/Fe] trend. However, we selected the stars from the synthetic survey to mimic the constraints in observations which we ultimately draw comparisons to.

Selection of stars for analysis
In this work, we aim to compare the age-[X/Fe] relations for the thin disk stars in the simulated Milky Way to those measured in the observed Milky Way lowα disk (e.g. Bedell et al. 2018, that have in-situ origins and have been the subject of most age-[X/Fe] studies. With this in mind, we have constrained our analysis of stars in the simulations to those that (1) have R<25 kpc, (2) [Fe/H] ≥ -0.25 dex, and (3) | Z |< 500 pc. We do this selection to obtain predominantly in-situ stars in the current day spatial thin disk, though we note that there is no explicit removal of ex-situ stars. We explain these criteria in detail below. We note that for ease of processing, we utilize the subsampled independent survey that is 1/100 the size of the full synthetic survey for m12i at LSR2.
The first criterion (1) restricts to stars within 25 kpc of the galactic center. This is motivated by Figure 3 of Sanderson et al. (2018), which shows the current locations of stars in the FIRE simulated galaxies colored by their formation locations. For m12i, these authors found that at R∼20 kpc, there is a sharp transition between stars that formed in-situ and stars that formed at greater distances, i.e. R >30 kpc. Azimuthally, this transition radius hovers around the R = 20 kpc value, though sometimes exceeds this range. We therefore opted to use R = 25 kpc as our radius cut to take into account the upper bound of this transition region. Interestingly, Sanderson et al. (2018) also found that this transition region between mostly in-situ vs mostly accreted stars is also traced by a transition in [Fe/H], which is a promising diagnostic for observational studies. This is also similar to the selection made in Bellardini et al. (2022) in defining the in-situ component of Milky Way-like galaxies in FIRE-2.
As a sanity check for this selection, we investigate how the spatial distributions, specifically in current (cylindrical) Radius, R, and current height from the disk plane, Z, of our resulting sample change as a function of age and metallicity. This is also motivated from observations Yuxi et al. 2021) that the tight relationship between [X/Fe] and age at fixed [Fe/H] is due to stars reflecting the chemistry of their location; therefore, we expect these distributions to be distinct .25 dex to -2.5 dex) as well, but focus on these three, which are sufficient to succinctly portray the differences in the R and Z distributions at high and low [Fe/H]. Each R and Z distribution is broken down into the separate contributions from different stellar ages, i.e., <1 Gyr (light blue), 1-3 Gyr (dark blue), 3-5 Gyr (green), 5-7 Gyr (olive), 7-10 Gyr (gold), and >10 Gyr (brick red). We also mark the median R and Z values with the purple dashed line. Though our first criterion makes a cut at R<25 kpc to select the predominantly in-situ population, not many stars lie in that region. In order to better compare the R trends at different ages, we therefore zoom in and only show the distributions up to R = 20 kpc.
At solar metallicity, the R distribution peaks at R = 8 kpc, slightly coinciding with the median R marked by the purple dashed line at R = 7.5 kpc. Meanwhile the R distributions for the lower metallicity samples (i.e. [Fe/H] = -0.25 dex and -0.75 dex) peak close to, but not quite at the center of the galaxy at R = 2 kpc, in line with the results from El-Badry et al. (2018) for the oldest stars in the FIRE-2 simulations. The median R for both low metallicity bins is higher than the peaks in their distributions and lies at R = 4 kpc, as both distributions have long tails towards larger R. For the bin that has a size of ±0.025 dex, we mark the distribution of stars at different ages i.e. <1 Gyr (light blue), 1-3 Gyr (dark blue), 3-5 Gyr (green), 5-7 Gyr (olive), 7-10 Gyr (gold), and >10 Gyr (brick red). We also mark the median R and Z at these metallicities (dashed purple line). The R distributions move from the inner to the outer regions while the Z distributions change from broad to narrow as we go from older to younger stars. Z distributions, the trend at solar metallicity resembles a normal distribution with mean at Z = 0. At lower metallicities, however, the distributions deviate from a Gaussian and exhibit two peaks, below and above the plane of the disk. The dispersion in Z also increases from high to low [Fe/H], as indicated on the upper left corner of each subplot in Figure 2, from σ Z = 1.10 kpc at [Fe/H] = 0, to σ Z = 2.05 kpc at [Fe/H] = −0.25 dex, and to σ Z = 3.27 kpc at [Fe/H] = −0.75 dex. This trend of increasing Z dispersion with increasing age and decreasing metallicity exists in the Milky Way as well (Meusinger et al. 1991), and could be explained by stars being born in a thicker disk at early times (Wisnioski et al. 2015;Ma et al. 2017;Bird et al. 2021;Yu et al. 2021), or by being heated kinematically through interactions with giant molecular clouds, spiral arms, and bar (Aumer et al. 2016) and with satellites (Hopkins et al. 2008;Villalobos & Helmi 2008Sales et al. 2009;Laporte et al. 2019). On the other hand, the double-peak in the Z distribution and the asymmetry about Z = 0 is due to the fact that it is harder to see stars through the disk because of the dust and that older stars have larger scale heights.
In H] = -0.25 dex bin is comprised of stars of different ages, the dominant stellar age is 5-7 Gyr (olive green line); therefore, the total R distribution (gray), across all ages, consequently follows the R distribution of the 5 − 7 Gyr stellar population. We note, however, that within the same [Fe/H] bin, these R (and Z) distributions at different ages highly overlap.
At [Fe/H] = 0, the R distribution peak at higher R values for younger stars than older stars. This is less evident for the lower [Fe/H] bins, where there are fewer younger stars. As seen for the [Fe/H] = −0.25 dex sample, the peak of the 5 − 7 Gyr distribution has R value very close to that of the 7 − 10 Gyr distribution. At the lowest [Fe/H] bin, the R where the distributions peak for stars with ages 5-7 Gyr, 7-10 Gyr, and >10 Gyr all coincide at ∼2 kpc, showing no sequence in R with the dominant ages. Therefore, we primarily observe evidence of an inside-out radial formation for younger and higher metallicity stars, i.e. ages < 7 Gyr and [Fe/H] > -0.25 dex. That is, only with these cuts do we isolate an inside-out forming disk with a clear current day age gradient at a fixed [Fe/H]. This age is broadly consistent with the transition age found by Yu et al. (2021) when FIRE-2 galaxies change from a bursty to smooth star formation mode. This is also the transition age that Bellardini et al. (2022) found for their sample of FIRE-2 Milky Way-like galaxies, when the stellar abundances changed from being dominated by azimuthal variations to radial variations. All Z distributions show increasing dispersion with increasing age at a given [Fe/H]. For the [Fe/H] = 0 dex sample, the Z distributions for stars with ages <1 Gyr and 1-3 Gyr are Gaussian, but the distributions for older stars deviate from this behavior. This is more obvious at [Fe/H] = -0.25 dex and [Fe/H] = −0.75 dex, which are dominated by older stars that show this double peak more, as the effect of not being able to "see" through the disk is more significantly felt. These R and Z trends that we find in the simulations are in congruence with the inside-out and upside-down scenario for the observed Milky Way stellar age-velocity relation as shown by Bird et al. (2021), and is further support for the Milky Waylike nature of this simulation and its earlier version in FIRE-1 (e.g., Ma et al. 2017).
With this exploration and these trends in mind, we therefore apply a second criterion in investigating the age-[X/Fe] relations (see Section 5) and (2) only select stars with [Fe/H] > -0.25 dex to ensure that we are probing the regime in metallicity in which chemical evolution has been temporally smooth, and age and [Fe/H] are a link to a star's location.
Lastly, (3) we spatially constrain our analysis to the simulation's thin disk, which is the focus of most age-[X/Fe] trend studies in observations (e.g. Nissen 2015, Bedell et al. 2018. We selected stars with Z within ±500 pc of the disk plane motivated by the measured scale height of the thin disk for m12i (Sanderson et al. 2020) and by the Z distributions of the radiallydistinct mono-age, mono-abundance populations shown in Figure 2.
As a caveat, we note that this comparison is not entirely consistent because in observational studies, the thin disk selection is usually done in the [α/Fe] vs [Fe/H] plane, while in most simulations and in this study, the criteria for the thin disk is in approximation a spatial selection. However, for the purpose of our analysis, this means we compare the stars in the disk with small vertical oscillations around the plane, in both simulation and data.

OBSERVATIONAL DATA
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, σ 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 < T eff < 7000 K. We further apply a cut in [α/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-α 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 σ 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 ≈ 15,000 red clump field stars from Ness et al. (2019) in the low-α disk, where the authors utilize the apogee sur- Kernel density estimations (KDE) for R birth of stars within age bins centered at 1 Gyr (purple) to 9 Gyr (yellow) in steps of 2 Gyr, prior to the [Fe/H] cut. The KDEs are normalized to the number of stars per age group. Stars of different ages trace distributions that, though highly overlapping, show distinct birth annuli. However, the 9 Gyr stellar population does not follow this trend due to a more chaotic mode of star formation at this point in the m12i's history.
vey'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 σ int as the measured scatter from the best fit line will be a combination of the σ int and the error in the abundance measurement, σ 2 abund . That is: σ 2 total = σ 2 abund + σ 2 int .

THE SPATIAL DISTRIBUTION OF STELLAR CHEMICAL ABUNDANCES AT BIRTH
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, 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 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 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 ∼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 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 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 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 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 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 ∼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 α-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-α disk taken from galah DR3 written in brackets. We note however that the age bins used for the observations are much bigger (e.g., ±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.  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 14 N, decrease in 12 C, and slight decrease in 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 = 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 Bensby et al. 2004), and loosely, age, as is expected of an α-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.
α-elements (Magnesium, Silicon, Sulfur, Calcium): Lastly, we consider the other α-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 α-elements, which are mainly produced through SNe II with shorter timescales, and diluted via SNIa that have delayed and longer timescales, therefore decreasing the α-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 α-elements. Mg, like O, is a hydrostatic α-element produced during the shell burning phase in massive stars. Si, S, and Ca are explosive α-elements that also have non-negligible contributions from SNIa. For the explosive α-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 α-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 R < 25kpc trace birth locations (Figure 4) that overlap quite broadly but peak at different R birth and exhibit [Fe/H] and [X/Fe] distributions that also show distinct trends (Figure 5 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) 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.

AGE-ABUNDANCE-PRESENT-DAY-LOCATION TRENDS AT FIXED [FE/H]
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, R now compared to 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 R now and 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.  . Current and birth radius for stars of different age ranges. Kernel density estimations for Rnow (solid) vs R birth (dashed) of stars within age bins centered at 1 Gyr (purple) to 9 Gyr (yellow) in steps of 2 Gyr (top to bottom). The KDEs are normalized to the number of stars per age group. The location KDEs, both for Rnow and R birth , exhibit similar tracks for stellar populations of the same age. This shows that although stars move from their birth location as shown by the changing R distributions, the bulk expectation for an inside-out galaxy growth is still observed, with old stars concentrated towards the center and young stars in the outskirts.  from these fits are included in each figure sub-panel i.e., the median intrinsic dispersion (σ int ), the slope, m, and the intercept, b. Similar to the age binning derivation in Section 4, we determined the best binning in [Fe/H] by progressively calculating the σ int at different bin values and choosing the bin that gives the smallest value so that its measurement is the least sensitive to the size of the [Fe/H] bin. We note that observational studies of the age-[X/Fe] relations typically have fixed and larger [Fe/H] bins (e.g., 0.05-0.10 dex) compared to this work, but we have added the step of determining the best binning to remove the dependence of the σ int in [X/Fe] on the scatter in [Fe/H]. Nonetheless, the discussion on the σ int from simulations and observations in Section 6 show that even with different [Fe/H] binning, their comparison is still viable. We excluded the running value of the intrinsic dispersion, σ int , calculated for age-[X/Fe] bins with fewer than 500 stars, to avoid being skewed to larger σ int due to sparsity in certain regions of the age-[X/Fe] parameter space. However, we note that sparsity still affects some of our analysis.
The cuts we applied (see Section 3) allow for a more direct comparison of the age-[X/Fe] trends between simulations and observations. However, there are some discrepancies in the scaling of the [X/Fe] between observations and simulations. Given this, we focus more on the similarities and differences in the [X/Fe] trends with age. Absolute differences in the mean element abundances between data and simulation are not unusual and could be caused by many things in the simulations, like the different star formation history in m12i compared to the Milky Way (Sanderson et al. 2020 Ness et al. 2019). As discussed in Section 4, C is produced mainly in SNe II in observations, and as a consequence, should have higher abundance for older stars than younger stars. In fact, Nissen (2015) found that [C/Fe] has a larger slope with age than α-elements, and proposed that a top-heavy IMF can possibly explain the strong trend in C abundance with age. There is a small intrinsic dispersion around the age-[C/Fe] trend, of σ int =0.019 dex.
The age-[N/Fe] trend shows a decreasing trend with older ages, similar to what is found in observations of RGB stars from apogee ). The measured slopes from the Milky Way observations and the simulations agree, with m=-0.019 dex Gyr −1 from this study (and σ int =0.019 dex) and m=-0.02 dex Gyr −1 from Ness et al. (2019). Unfortunately, N is hard to measure for stars in the optical so we are unable to compare to more observational studies.
The [O/Fe] trend with age is slightly decreasing, but essentially flat, having m=-0.001 dex Gyr −1 and a σ int =0.014 dex around the best fit line. This is contrary to what is observed for the Milky Way (Nissen 2015;Bedell et al. 2018;Ness et al. 2019;Sharma et al. 2022). O is a hydrostatic α-element dispersed mostly through SNe II and is expected to decrease with younger stars. As noted in Section 4, the metallicity-dependence of O yields from stellar winds in the simulations significantly contribute at higher [Fe/H], therefore driving the trend to be flat or even decreasing with older age.
All the remaining α-elements-Mg, Si, S, and Cashow positive age-[X/Fe] trends, in line with expectations from Galactic chemical evolution work (Kobayashi et al. 2020) and from observations (Nissen 2015;Bedell et al. 2018;Ness et al. 2019;Sharma et al. 2022) with the exception of Ca, where observations find a flat age-[Ca/Fe] trend. Ca is an explosive α-element, which is expected to have a similar slope as compared to other explosive α-elements like Si and S. We observe this in the simulations where the slopes for Si, S, and Ca are 0.006, 0.006, and 0.005 dex Gyr −1 respectively, as the SNIa contributions flatten the trend with age, while in comparison, the hydrostatic α-element Mg has m=0.008 dex Gyr −1 , a steeper trend because it is mostly produced in SNe II. Ca, aside from showing a flat relationship with age in observations also shows a flat relationship in stellar velocity dispersion with increasing abundance, contrary to what is found for other α-elements (Conroy et al. 2014). This is especially intriguing because there is a known age-velocity dispersion relation in observations wherein stellar populations with older ages have higher stellar velocity dispersion. One way to explain this curious trend of Ca with age is to include contributions from other sources, such as a subclass of supernovae called "calcium-rich gap transients" that would flatten the trend at later times (Mulchaey et al. 2014;Nissen 2015).
For all the elements except Fe, the running σ int gets larger at older ages where there are (1) fewer stars in general at solar metallicity and (2) the stars are born from less homogeneous gas.

Solar neighborhood at [Fe/H]= -0.25 dex
We apply the same spatial selection and binning determination as per  Figure 8 show more clumpiness and scatter than the [Fe/H] = 0 dex trend, especially at older ages i.e., one clump at 5<age<7.5 Gyr and another at age >8 Gyr. For each element, the age-[X/Fe] trend is steeper compared to the [Fe/H] = 0 dex sample (see Figure 9 for comparison), even without the inclusion of stars in the oldest clump. In fact, the age-[O/Fe] trend changes from slightly decreasing (or flat) at [Fe/H] = 0 dex to slightly increasing with a slope of m=0.002 dex Gyr −1 , due to the less significant contribution from stellar winds at this point of m12i's chemical enrichment. The median σ int values are all higher for all the elements in this lowmetallicity regime, except for Fe, where we made a cut in the sample selection. Ness et al. (2019)  In this section, we explore the age-[X/Fe] trends in the simulations for stars at [Fe/H] = 0 ± 0.025 dex, Z <500 pc, and at different distance slices to represent the inner disk (1-3 kpc), outer disk (13-15 kpc), and the whole disk (<25 kpc) as shown in Figure 10. We particularly want to understand if different locations in the galaxy exhibit distinct age-[X/Fe] trends, and how this ties to the m12i's history.  Compared to the solar neighborhood trend in Figure  7, the inner disk region (top set of subplots) shows higher σ int around the age-[X/Fe] relations. There is also a more pronounced increase in σ int for ages > 4 Gyr. In fact, this behavior extends to much older ages where the number of stars is below the threshold that we have set. These old stars, as we have shown in Section 4, and Figures 4 and 5, are born with a larger spread in their abundances when the star formation in the galaxy was more chaotic and clumpier. The stellar populations and the general age-[X/Fe] relations in the inner disk are both more complex, a reflection of the inner disk region's intricate history where stars of different ages all contribute. Nonetheless, we adopt a linear fit to the age-[X/Fe] distributions so as to enable more direct and interpretable comparison to observations and to the age-[X/Fe] trends at other locations. The increase that we see in σ int , which signals the difference between steady and bursty star formation, is supported by the recent study by Yu et al. (2021), which finds a similar transition at ∼3 Gyr for m12i. This is only seen closer to the inner disk, but is slightly noticeable in the solar neighborhood sample as well.
At the same [Fe/H] (i.e. [Fe/H] = 0 dex), the inner disk shows age-[X/Fe] relations in the α-elements that are higher at older ages and lower at younger ages compared to the solar neighborhood ( Figure 7) and outer disk (middle set of subplots in Figure 10). The higher α-element abundance in the inner disk is in agreement with the rapid star formation and enrichment scenario in the center of galaxies, before the onset of SNIa. On the other hand, the lower α-element abundance in the same region is in line with the longer history of star formation in the inner regions of the galaxy where the [α/Fe] abundances have been more diluted by SNIa compared to the outskirts of the galaxy. Compared to the age-  Figure 7. The top set of plots are stars in the inner disk at 1< R <3 kpc, the center set are for stars in the outer disk at 11< R <13 kpc, and the bottom set are for all stars within R = 25 kpc and Z±500 pc. The scatter around the age-[X/Fe] relations is higher for the inner disk compared to the outer disk.This is a consequence of older stars (of which there are more in the inner disk) exhibit wider abundance distributions. Interestingly, when all distance slices are included, the median σint is still quite small at σint<0.025 dex.
[X/Fe] trends in the inner disk, the stars in the outer disk are mostly young with ages <5 Gyr and therefore have a more well-behaved age-[X/Fe] trend across the board that is well-fit by a linear trend.
Lastly, we look at the age-[X/Fe] trends for the entire disk (bottom set of subplots in Figure 10). Even with the inclusion of stars from various locations in the galaxy-from the inner disk that experienced different modes of star formation to the outer disk that has expe-rienced later and more steady star formation-the median σ int remains small with values of σ int <0.025 dex. We note that the gray data points are meant to show the total extent of the age-[X/Fe] trends, but not the density. This could be a misleading representation for the whole disk, where at ages <5 Gyr, the running σ int is small, but the extent in [X/Fe] at a given age is large. Nonetheless, this provides a direct comparison to the age-[X/Fe] relations discussed earlier. We see an in-crease in the running σ int towards older stars, bearing similarity to the inner disk. The small σ int measured in m12i, even after including all the stars "observed" across the full disk sample points to a case where this small σ int is representative of a population whereby knowing the age and [Fe/H] of a star can predict its abundance, with increasing scatter at older ages. Ultimately, the ability to predict the abundance from the age and [Fe/H], highly depends on when the galaxy changed to a steady star formation from a more bursty mode as this determines the σ int and precision.

INTRINSIC DISPERSION IN THE AGE-[X/FE] -METALLICITY RELATION ACROSS THE GALAXY'S DISK
We take a deeper look at the σ 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.

σ int in simulation vs observations
The intrinsic scatter around the age-[X/Fe] trends, σ int , as measured for stars in the solar neighborhood at solar [Fe/H], from this work and other observational studies, is shown in Figure 11. The diamonds are stars from the simulations (Ananke), while the circles are solar twins (Bedell et al. 2018), and crosses are RGB stars ) from observations. As noted in Section 3, these observational data take into account the effects of the [X/Fe] and age uncertainties on the resulting σ int . The symbols are colored by the age-[X/Fe] slope where red is negative and decreasing with age and blue is positive and increasing with age.
The σ int for the elements produced in the simulated Milky Way are similar to the σ int from Milky Way observations, both for the solar twins and RGB stars, and are generally σ int ≤0.02 dex, marked by the gray dashed line. This is especially true for the α-elements Mg and Si. This small scatter is similarly found by Bellardini et al. (2022) for [Mg/Fe] at different age bins in their sample of FIRE-2 Milky Way-like galaxies. The general similarity in the σ int means that the level of homogeneity of star forming gas in the Milky Way at a given metallicity is well-reproduced in the simulated galaxy m12i.
As discussed in Sections 4 and 5, C, O, and Ca exhibit dissimilar trends with age compared to observations (and again note we refer to age-[X/Fe] relations in narrow bins of [Fe/H]). The element C shows a negative slope with age in the simulation, while it has a positive slope measured both from dwarf and giant stars in the Milky Way. The element Ca, on the other hand, has a positive slope, typical of α-elements in the simulations. However, Ca has an observed negative slope with age. The element O is essentially flat with age in the simulations, in contrast to the increasing trend with age seen in observations. In a simple exploration of the age-[X/Fe] relations in the other galaxies in Ananke (m12m and m12f) with different star formation histories, we find that they similarly exhibit these discrepancies in the direction of the trends. This highlights the need for updates to the chemical enrichment sources, yields, and rates prescribed in the star formation in cosmological simulations in order to further their utility for Galactic archaeology purposes.
There have been many recent efforts to look at the age-[X/Fe] trends in different Milky Way stellar populations or spectral types. These recent efforts are a reflection of how powerful a tool chemical abundances could be, especially for deriving ages (e.g. Hayden et al. 2020). Yuxi et al. (2021) compared the age-[X/Fe] trends for the high-and low-α disk using apogee red clump stars and found that the median σ int for both stellar populations are similar at σ int ∼0.04 dex; even with stellar populations with different chemical enrichment histories, the abundance scatter is comparable. Although there is no α-element dichotomy in m12i in the [α/Fe]-[Fe/H] plane, we do find that the σ int measured at different radial parts of the disk are small and generally σ int <0.04 dex, even if they experienced different formation histories. In this sense, our work agrees with the findings from Yuxi et al. (2021). However, a further exploration of the age-[X/Fe] relations for different galaxy components, as Nikakhtar et al. (2021) recently determined for the Ananke galaxies, would be worth exploring in detail to compare to observations. Using galah DR3 data, Sharma et al. (2022) found that for ages <10 Gyr, a disk star's abundance can be predicted from its age and [Fe/H] with a σ int of 0.03 dex. They also found that comparing the slopes in the age-[X/Fe] trends to those from the [Fe/H]-[X/Fe] relations split elements into three groups, corresponding to their main channel of production: SNe II, SNIa, and stellar winds; these are the same sources used for the chemical enrichment in our simulated Milky Way.

σ int at different disk locations and [Fe/H] in the simulations
We now explore the σ int of each element around their age-[X/Fe] relations for the seven elements included in this work (excluding [Fe/H]) at different radial locations in the disk (e.g. 1-3 kpc, 4-6 kpc, 7-9 kpc, and 10-12 kpc) and at different metallicities (solar at [Fe/H] = 0 dex and low at [Fe/H] = -0.25 dex) as shown in Figure  12. We color our markers by the median age of that

stellar population in a given location and [Fe/H] bin.
The 2kpc radial bin is consistent with the radial range over which the stars are similar in chemistry (Bellardini et al. 2021 vary with location, and at the same radius bin, the difference in median ages is not proportional to the difference in σ int (for example, comparing between the 5 kpc and 11 kpc bin). Ultimately, there is the caveat that the σ int trends as a function of [Fe/H] and location in this work are unique to m12i. 6.3. σ int in the real solar neighborhood at different [Fe/H] In addition to the results from the simulation data, we derive the scatter around age-[X/Fe] relations, σ int , from galah DR3 observations of real solar neighborhood stars at different metallicities, shown in Figure 12. The following cuts were made to ensure that the stars from galah data (described in Section 3) are comparable to our bins of simulated stars and that the abundances are of good quality: • 7<R<9 kpc, • signal-to-noise ratio >40, • flag sp=0, • flag X fe=0, • 5000<T eff <6100 K, and where X = C, O, Mg, Si, and Ca-the elements this study has in common with galah. Note that our T eff and log g cuts are adopted from the MSTO stars cut from Sharma et al. (2022). With these selections, we end up with 2,882 stars for the solar [Fe/H] sample and 1,546 stars for the low [Fe/H] sample. In general, we performed a similar exercise as was done for the simulations, and calculated the median running dispersion in the best fit linear age-[X/Fe] trend for these elements.
Like the other comparison observational data (see Section 3), the uncertainties in [X/Fe] and age have to be taken into account to calculate the σ int . To account for the abundance uncertainties, we subtracted the mean σ 2 abund from the measured σ 2 total to get the σ 2 int . The typical value for σ abund is 0.05 dex for Mg, Si, and Ca and 0.10 dex for C and O. Next, to understand how the age uncertainty affects the measured σ int , we drew a new age from a normal distribution centered on the age given by BSTEP with the age uncertainty as standard deviation. The mean age error for our sample is 0.7 Gyr. We do this 100 times and find that the change in σ int is 0.01 dex and therefore negligible. to what we find in the simulations. All measured σ int from the observations are lower than 0.07 dex, and the difference in σ int between the two [Fe/H] bins ranges from <0.01 dex (e.g. C, Ca) to ∼0.03 dex (e.g. Mg). The σ int trend for Si in observations resembles the σ int in the simulations the most, while the σ int for Ca in galah resembles the simulations the least, where the σ int from the observations are larger by ∼0.03 dex compared to their simulated counterpart. Even with these differences, it is noteworthy that the σ int from the observations and simulations are very similar and of the same magnitude. This comparison highlights that even if the absolute values and trends for [X/Fe] in the simulations are dependent on many assumptions (i.e., on yields, rates, sources) and the star formation history, the general abundance scatter in observations is reproduced well (Hopkins et al. 2018;Escala et al. 2018;Bellardini et al. 2021Bellardini et al. , 2022. We expect that this agreement holds true even if realistic observational uncertainties were modeled for m12i as (1) the effect of the abundance error is removed when deconvolving the sources of the scatter and (2) the slopes in the age-[X/Fe] relations in the simulations are shallower than what the observations measure. Applying age uncertainties in the simulations will therefore have negligible effects on the resulting σ int .

DISCUSSION AND SUMMARY
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 ∼0.02 dex. Indeed, this has been studied in greater detail and in many axes e.g., high vs low-α disk, high vs low [Fe/H], MSTO vs RGB stars, thanks to large spectroscopic survey data 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, metallic-ity, 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, σ int , generally reflect m12i's star formation history. Focusing on the solar [Fe/H] sample, at all distances, we find that the running σ int across age is constant at young ages, up to about ∼4 Gyr, and then increases to higher σ 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 ∼3 Gyr. Indeed, for stars that are <3 Gyr, the σ int is constant, for 3<age<6 Gyr, the σ int increases, and at >6 Gyr, the σ 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 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∼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 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 signif-icant overlap in their distributions of [X/Fe], 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 σ 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 σ int around the age-[X/Fe] relations from this work to those derived in observations. Specifically, in Figure 11 we show the σ int in the solar neighborhood at solar metallicity for the stars in Ananke, solar twins (Bedell et al. 2018), and RGB stars . We find that the σ int for m12i and the Milky Way are similar, and in fact, some elements such as Mg and Si have near-identical σ 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 α 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 α-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 α 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 Waylike 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 σ int for the different Milky Way galaxies in Ananke that have masses 1-2 ×10 12 M but have varying morphologies and star formation histories. This would be important in contextualizing what brings about the observed σ int in the Milky Way by comparing to the derived σ int in simulated Milky Way galaxies with different properties. In line with this, we have shown in Section 6 how the σ 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 σ 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 σ 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.

ACKNOWLEDGMENTS
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 sup-port 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.

A. YIELD PRESCRIPTIONS IN THE SIMULATION
As discussed in Section 2, the stellar evolution in FIRE-2 Hopkins et al. (2008) is produced from STARBURST99 (Leitherer et al. 1999) with a Kroupa (2001) IMF. The three sources of nucleosynthesis included are SNe II, SNIa, and stellar winds from OB/AGB stars. In the simulations, SNe II rates are a function of time (i.e. age of a star particle) dictated by the following:  Table A1. Lastly, mass-loss from OB/AGB stellar winds were included whose yields are listed on the third column of Table A1. For more details on the stellar evolution in FIRE-2, we refer to Appendix A in Hopkins et al. (2018). Table A1. Element yields from different sources. SNe II yields are from Nomoto et al. (2006), SNIa from Iwamoto et al. (1999), and stellar winds from the compilation by Wiersma et al. (2009) of van den Hoek & Groenewegen (1997, Marigo (2001), and Izzard et al. (2004). Yields are in terms of M . We note that the N abundance has a metallicity dependence and is multiplied by N where N =max(Z/Z ,1.65). In Section 4, we determined the optimal binning in age to look at the R birth , current R, [Fe/H] and [X/Fe] distributions as a function of age. We started with a bin size of ± 0.5 Gyr and scrutinized the [Fe/H] vs R birth of the stars encompassed by this selection. We reduce the bin size by half, and calculate how the dispersion across the linear fit changes. A decrease in dispersion warrants a further decrease in the age bin; an increase in dispersion or a sample size <50,000 stars stops the iteration. In general, we define the age bin as sufficient when it is (1) not too wide such that the measured dispersion in the [Fe/H] vs R birth trend is independent of the binning and (2) not too narrow to be affected by small number statistics. Figure B1 illustrates this process for 5 Gyr old stars, wherein we show their 2D histogram, with the best fit line for [Fe/H] vs R birth , and the running dispersion around this trend. We also note the number of stars included in the bin in the bottom left corner and the median running dispersion in the top right corner.

Elements
As we reduced the age bin size from 0.5 Gyr, the dispersion in [Fe/H] vs R birth also continue to decrease until we reach 0.03 Gyr. At this point, the nth iteration has a small sample size that drives up the dispersion (e.g. 0.147 dex). We therefore choose the age binning from the n-1 iteration which gives us the smallest dispersion in the [Fe/H] vs R birth diagram. We perform a similar exercise for the other ages included in the analysis.