© The Authors 2024

1. Introduction

The ultimate goal of Galactic archaeology – the study of the history of star formation as recorded in stellar properties – is to understand the formation and evolution of the Milky Way. To do so, though, we have to interpret its history from a present-day snapshot, which includes stellar kinematics, ages, and chemical abundances. Knowledge of the stars’ time and place of birth would allow us to directly trace the evolution of the Galaxy over time.

It is well established that stars radially migrate away from their birth sites (e.g. Sellwood & Binney 2002; Schönrich & Binney 2009; Minchev & Famaey 2010; Frankel et al. 2018; Carr et al. 2022), meaning that a given location in the Galaxy can be comprised of stars born at a variety of different radii depending on stellar age and migration rate (Minchev et al. 2018; Agertz et al. 2021; Carrillo et al. 2023). It has also been shown that radial migration affects stars relatively early in their lifetime (Kubryk et al. 2013; Grand & Kawata 2016; Feltzing et al. 2020), and therefore this mixing of stars from different birth radii (R_birth) at a given current location makes it challenging to unravel evolutionary events; radial abundance gradients become weaker (Pilkington et al. 2012; Minchev et al. 2012a, 2013; Kubryk et al. 2013; Ratcliffe et al. 2023) and chemo-kinematic relations can become blurred (e.g., Minchev et al. 2014).

Due to gravitational interactions with spiral arms, the central bar, and infalling satellites, stellar angular momentum can be permanently changed, making it impossible to trace a star’s position back to its birth site by simply using its current kinematics (Sellwood & Binney 2002; Minchev & Quillen 2006; Roškar et al. 2008; Quillen et al. 2009; Minchev & Famaey 2010; Minchev et al. 2012b; Di Matteo et al. 2013; Khoperskov et al. 2020). Chemical abundances, on the other hand, are an archaeological record of the physical conditions in the interstellar medium at the birth time and position of stars as they are relatively unchanging with time (Tinsley 1979). The idea of strong chemical tagging utilizes these observable abundances to recover dispersed star clusters, which are believed to be chemically homogeneous (Freeman & Bland-Hawthorn 2002; Bland-Hawthorn et al. 2010). While recent works have shown that strong chemical tagging is infeasible with current-day observational uncertainties and sample size (Lindegren & Feltzing 2013; Ting et al. 2015; Casamiquela et al. 2021), a weaker version of chemical tagging to recover stars born at similar times and places in the Galactic disk is more likely to succeed (Ratcliffe et al. 2020, 2022).

Recent works leverage chemical abundances of disk stars to estimate the effect of radial migration in the Milky Way disk (e.g. Minchev et al. 2018; Frankel et al. 2018, 2019, 2020; Lu et al. 2022a; Ratcliffe et al. 2023, 2024; Wang et al. 2024). One way to understand this effect and recover the Milky Way’s evolution is by estimating stellar birth radii from a star’s metallicity and age. These R_birth methods rely on a linear relationship between [Fe/H] and R_birth at a given look-back time, where the slope (birth metallicity gradient over cosmic time; ∇[Fe/H](τ)) and intercept (metallicity at a specific radius R over time; [Fe/H](R, τ)) are unknown. Through the development of a largely model independent approach, Minchev et al. (2018) simultaneously recovers stellar R_birth, ∇[Fe/H](τ), and [Fe/H](R = R_⊙, τ) by requiring that R_birth distributions of mono-age populations remain non-negative and follow expectations of inside-out disk formation (Matteucci & Francois 1989). With their discovery that the scatter in metallicity across age bins is linearly correlated with ∇[Fe/H](τ), Lu et al. (2022a) allow for the determination of R_birth to be predominantly self-contained, with the largest modeling requirement being that the evolution of [Fe/H] in the Galactic center smoothly and continuously enriched over time.

In recovering ∇[Fe/H](τ), Lu et al. (2022a) finds that the radial [Fe/H] gradient had a steepening event 8 − 11 Gyr ago, which they attribute to the Gaia-Sausage-Enceladus (GSE; Belokurov et al. 2018; Helmi et al. 2018; Haywood et al. 2018. See Prantzos et al. 2023 for an alternative argument). Applying the same method to higher spectral resolution data, Ratcliffe et al. (2023) find two additional steepening events about 4 and 6 Gyr ago, which coincided with passages of the Sagittarius dwarf galaxy (Ibata et al. 1994; Law & Majewski 2010). The effect of mergers on the radial metallicity gradient is investigated using cosmological hydrodynamical simulations in Buck et al. (2023), where it is found mergers brought in pristine gas and lowered the metallicity in the outer disk, therefore steepening the gradient (though potentially not over long time scales; Renaud et al. 2024). It is also shown that if the merger was large enough, the metallicity at the galactic center would decrease, in direct conflict with the monotonic [Fe/H] enrichment at the Galactic center used in Milky Way disk R_birth derivations (e.g. Buck 2020; Renaud et al. 2024).

Modern hydrodynamical cosmological simulations – such as NIHAO-UHD (Buck et al. 2020), Feedback In Realistic Environments-2 (FIRE-2; Garrison-Kimmel et al. 2018), VINTERGATAN (Agertz et al. 2021; Renaud et al. 2021a,b), and HESTIA (Libeskind et al. 2020; Khoperskov et al. 2023a) – contain a wealth of information for investigating the effect of different evolutionary histories. In particular, TNG50 (Pillepich et al. 2019; Nelson et al. 2019a,b) contains a sample of 198 Milky Way and Andromeda-like galaxies that have diverse galactic properties and histories (Pillepich et al. 2023). With access to a full set of stellar properties – including those that are not measurable for Milky Way stars (e.g. R_birth) – for a large variety of simulated Milky Way-like galaxies, we can investigate the dependencies of different evolutionary histories on our ability to recover stellar birth radii.

Previous works have only tested the ability to derive birth radii on a handful of simulated galaxies, with little to no investigation into how the central metallicity evolves with time. In this work, we further tested the assumptions required to estimate birth radii using the large Milky Way and Andromeda analogue sample from TNG50. Thanks to the large quantity of galaxies in this sample, we were able to test the validity of the assumptions on a diverse set of formation and evolutionary histories and explore which conditions are potentially violated in the Milky Way. This paper is organized as follows. Discussions of the data and R_birth method used for this work are given in Sections 2 and 3. Sections 4 and 5 present our findings on the evolution of the metallicity gradient and central metallicity assumptions, respectively. The discussion of our work in context with the Milky Way is given in Section 6. Section 7 gives our conclusions.

2. Milky Way and Andromeda analogues from TNG50

The TNG50 cosmological simulation (Pillepich et al. 2019; Nelson et al. 2019a,b) is the highest resolution run of the IllustrisTNG project (Pillepich et al. 2018; Marinacci et al. 2018; Naiman et al. 2018; Springel et al. 2018; Nelson et al. 2018, 2019a), and was run with the moving-mesh code AREPO (Springel 2010). The simulation evolved gas, stars, dark matter, black holes, and magnetic fields inside a box of width 51.7 comoving Mpc from redshift ≈127 to redshift 0, with a resolution of m_baryon = 8.5 × 10⁴ M_⊙, m_DM = 4.5 × 10⁵ M_⊙, and a spatial resolution of the star-forming ISM gas of ∼100 − 140 pc (Pillepich et al. 2019). Chemical enrichment of the interstellar medium accounted for heavy element nucleosynthesis from core-collapse supernovae (SNII), Type-Ia supernovae (SNIa), and asymptotic giant branch stars (Pillepich et al. 2018; Naiman et al. 2018). The yield tables were taken from Karakas (2010), Doherty et al. (2014), Fishlock et al. (2014) (AGB), Portinari et al. (1998), Kobayashi et al. (2006) (SNII), and Nomoto et al. (1997) (SNIa) (see Table 2 in Pillepich et al. 2018 for more information). The mixing processes were treated naturally in AREPO, with lower rates of numerical diffusion than fixed grid codes (Springel 2010; Naiman et al. 2018).

Within TNG50 there are 198 galaxies that are similar to the Milky Way or Andromeda, and were selected using the following selection criteria at redshift 0 given in Pillepich et al. (2023):

• The stellar mass of the galaxy is in the range M_*(< 30 kpc) = 10^{10.5 − 11.2} M_⊙;

• Disk-like morphology with spiral arms, chosen such that either the minor-to-major axis ratio of the galaxy’s stellar mass distribution is smaller than 0.45 or the galaxy appears disk-like and exhibits spiral arms through visual inspection of three-band images in edge-on and face-on projections;

• No other galaxy with stellar mass ≥10^10.5 M_⊙ is within 500 kpc and the total mass of the halo host is smaller than M_200c(host) < 10¹³ M_⊙.

We partner the Milky Way and Andromeda-like Satellites catalog¹ (Pillepich et al. 2023) with bar strength measurements from Khoperskov et al. (2023b) (taken as the peak value of m = 2 Fourier harmonics of the density distribution) and disk scale lengths from Sotillo-Ramos et al. (2022) (measured by fitting an exponential profile to the radial stellar surface density distribution excluding the bulge region) at redshift 0. Merger history information for each subhalo was recovered using SUBLINKGAL (Rodriguez-Gomez et al. 2015). The gas mass ratio of a merger is defined as the ratio of the gas mass of the secondary progenitor to the first progenitor at the time when the secondary progenitor achieves its maximum gas mass².

To make the best comparisons of this sample with the Milky Way, we scaled each galaxy to have the same disk scale length (3.5 kpc; Bland-Hawthorn & Gerhard 2016) and rotational velocity at the sun (238 km/s; Bland-Hawthorn & Gerhard 2016) as the Milky Way disk. We also performed the following selection cuts to capture the stellar disk: R_birth < 15 kpc, R < 15 kpc, |z|< 1 kpc, |z_birth|< 1 kpc, ecc < 0.5. Eccentricity was calculated using the following equation:

$$\begin{array}{c}\hfill \mathrm{ecc}=\sqrt{v_r^2+2{(v_{\varphi }-v_0)}^2}/(\sqrt{2}{v}_0),\end{array}$$

where v_r and v_ϕ are the radial and azimuthal velocities of the stellar particle, and v₀ is the rotational velocity of the galaxy at 2.2 disk scale lengths. We only used stellar particles that were born bound and are currently bound to the host galaxy to ensure stellar particles formed in situ. We also scaled the [Fe/H] abundance to the solar value from Asplund et al. (2009). Due to uncertainties in chemical yields and modeling in simulations (Buck et al. 2021), we focus our analysis on the relative values of [Fe/H] over time, and not the absolute values.

3. Birth radii Method

The goal of this paper is to extensively investigate the assumptions required to recover the evolution of the metallicity profile in galactic disks and additionally stellar birth radii. The recent method proposed by Lu et al. (2022a) (based on Minchev et al. 2018) assumes that for any lookback time (τ) we can write [Fe/H] as a function of the metallicity gradient at that time (∇[Fe/H](τ)), birth radius, and the metallicity at the Galactic center ([Fe/H](0, τ)):

$$\begin{array}{c}\hfill [\mathrm{Fe}/\mathrm{H}](R_{\mathrm{birth}},\tau )=\mathrm{\nabla }[\mathrm{Fe}/\mathrm{H}](\tau )R_{\mathrm{birth}}+[\mathrm{Fe}/\mathrm{H}](0,\tau ).\end{array}$$(1)

Then, by rearranging Equation (1), R_birth can be found as a function of metallicity and age:

$$\begin{array}{c}\hfill R_{\mathrm{birth}}(\mathrm{age},[\mathrm{Fe}/\mathrm{H}])=\frac{[\mathrm{Fe}/\mathrm{H}]-[\mathrm{Fe}/\mathrm{H}](0,\tau )}{\mathrm{\nabla }[\mathrm{Fe}/\mathrm{H}](\tau )}·\end{array}$$(2)

The main assumptions required can be broken down into two categories: assumptions regarding (1) the metallicity gradient and (2) the central metallicity. Our first results section (Section 4) investigates the assumptions surrounding the metallicity gradient, while our second results section (Section 5) explores the assumptions involved in the evolution of [Fe/H](0, τ). Section 6.1 puts this method in context with the Milky Way.

4. Results I: recovering the evolution of the [Fe/H] radial gradient

We start by examining the assumptions required to recover the evolution of the radial metallicity gradient with time. Namely, we tested the validity of a linear relationship between [Fe/H] and birth radius (Section 4.1) and uncovered which galaxies show a correlation between the metallicity gradient and range in metallicity ($\mathrm{Range}\stackrel{\sim }{[\mathrm{Fe}/\mathrm{H}]}(\mathrm{age})$) as first proposed in Lu et al. (2022a) (Section 4.2).

4.1. Linear relationship between [Fe/H] and birth radius

A necessary requirement in the ability to recover stellar birth radii is a linear birth metallicity gradient throughout the galactic disk. This is supported in observations of young stars in the Milky Way (Deharveng et al. 2000; Esteban et al. 2017; Arellano-Córdova et al. 2021) and is shown in some simulations (Vincenzo & Kobayashi 2018; Lu et al. 2022b). However, some simulations exhibit a non-linear relationship between [Fe/H] and R_birth, with the inner few kpc having a steeper slope (e.g. Agertz et al. 2021; Bellardini et al. 2022). Therefore, we begin our analysis by investigating the likelihood of a galaxy hosting a linear birth metallicity gradient in the TNG50 Milky Way and Andromeda analogues sample.

Figure 1 shows the median correlation (with the shaded region representing ±1 standard deviation) between [Fe/H] and R_birth over lookback time colored by bar strength as defined in Khoperskov et al. (2023b). Since the purpose of this work is to investigate our ability of recovering birth radii with the method used in Lu et al. (2022a), Ratcliffe et al. (2023, 2024), we chose to examine the correlation in age bins of 1 Gyr (similar trends are shown for age bins of 0.5 Gyr). To accommodate for outliers and strong migration during this period (Section 6.1.2), we looked at the correlation of the [Fe/H] and R_birth running median. Therefore, the correlation is defined as

Fig. 1. Median absolute correlation between [Fe/H] and Rbirth for galaxies in the TNG50 Milky Way and Andromeda catalog across lookback time grouped by bar strength, with the shaded areas showing ±1 standard deviation about the median. The correlation between [Fe/H] and Rbirth was calculated for age bins of 1 Gyr. To accommodate for outliers and higher feedback, the correlation was calculated for the running median, using Rbirth bins of 0.25 kpc. We do not find that bar strength has an impact on the linearity of the average birth metallicity gradient. The large variability in the correlation between [Fe/H] and Rbirth for stronger barred galaxies at smaller lookback times is predominantly due to weak metallicity gradients, even though the gradient itself is linear.

$$\begin{array}{c}\hfill \frac{{\sum }_{i=1}^n(R_{\mathrm{birth},i}-\overline{R_{\mathrm{birth}}})({[\mathrm{Fe}/\mathrm{H}]}_i-\overline{[\mathrm{Fe}/\mathrm{H}]})}{{\sum }_{i=1}^n{(R_{\mathrm{birth},i}-\overline{R_{\mathrm{birth}}})}^2{\displaystyle {\sum }_{i=1}^n}{({[\mathrm{Fe}/\mathrm{H}]}_i-\overline{[\mathrm{Fe}/\mathrm{H}]})}^2},\end{array}$$

where n is the number of radial bins used of width 0.25 kpc, [Fe/H]_i and R_birth, i represent the median [Fe/H] and average R_birth of the radial bin, and the bar denotes the average R_birth and [Fe/H] of that bin. Overall, we find a high correlation between [Fe/H] and R_birth regardless of bar strength, in agreement with previous works on linear birth metallicity gradients (e.g. Lu et al. 2022b). However, unlike Lu et al. (2022b), we find that a strong linear metallicity gradient exists already at early times when the photo-galactic disc is just starting to form.

While the median [Fe/H] gradient is linear across the TNG50 Milky Way and Andromeda-like galaxies of our sample, there are a number of galaxies which have a lower correlation between metallicity and birth radius as shown by the shaded regions in Figure 1. Many of these galaxies are given an artificially low correlation value due to their seemingly flat [Fe/H] gradients in the past few Gyr, even though their gradient is technically linear. Another key factor in our determination of linearity of the metallicity gradient is the number of disk particles. About 20 galaxies have few stellar disk particles with younger ages, with some having nearly no stellar particles with ages < 4 Gyr. This may be related to late massive mergers destroying the disk, which is found in 18 TNG50 analogues (Sotillo-Ramos et al. 2022).

Since a non-linear metallicity gradient prevents the recovery of birth radii, we only examined galaxies that have a correlation between R_birth and the [Fe/H] running median <  − 0.85. This corresponds to times when ∼75% of the variation in [Fe/H] is predictable from R_birth. To accommodate for weak gradients that are given low correlation values despite being linear, we choose to ignore those years when deciding how many recent Gyr the disk had a linear gradient. We, therefore, removed 17 galaxies (8.5% of our sample) due to their weak correlation between [Fe/H] and R_birth in the most recent Gyr for the remainder of this section. These 17 galaxies span a variety of bar strengths, and appear to be non-Milky Way-like in a variety of ways (see Appendix A for more information).

4.2. Correlation between the metallicity gradient and scatter

Lu et al. (2022a) shows that in the NIHAO-UHD (Buck et al. 2020) and Hestia simulations (Libeskind et al. 2020; Khoperskov et al. 2023a), there is a correlation between the scatter in [Fe/H] across different age bins and the evolution of the metallicity gradient with time. An example illustrating this correlation for a TNG50 galaxy is given in Figure 2. The density of the [Fe/H]–R_birth plane for stars born 2 − 3 Gyr ago for this example galaxy is depicted in the left panel, with the fitted gradient overlaid in black. The red dashed lines correspond to the boundary of the measured range in [Fe/H] during this age bin. The right panel shows ∇[Fe/H](τ) vs $\mathrm{Range}\stackrel{\sim }{[\mathrm{Fe}/\mathrm{H}]}(\mathrm{age})$ for the same galaxy, with each point corresponding to a different age bin. $\mathrm{Range}\stackrel{\sim }{[\mathrm{Fe}/\mathrm{H}]}(\mathrm{age})$ here refers to the normalized Range in [Fe/H] over time, where the time of the largest range corresponds to $\mathrm{Range}\stackrel{\sim }{[\mathrm{Fe}/\mathrm{H}]}(\mathrm{age})$ = 1 and the time of the smallest range has $\mathrm{Range}\stackrel{\sim }{[\mathrm{Fe}/\mathrm{H}]}(\mathrm{age})$ = 0.

Fig. 2. Example (subhalo 577372) illustrating the recovery of ∇[Fe/H](τ) from $\mathrm{Range}\stackrel{\sim }{[\mathrm{Fe}/\mathrm{H}]}(\mathrm{age})$. Left: Density distribution of the [Fe/H]–Rbirth plane for stars born 2 − 3 Gyr ago. The correlation between the [Fe/H] running median and Rbirth is given in the top right corner. The fitted true birth gradient is shown as the black line, and the maximum and minimum [Fe/H] used to determine the range are shown as the dashed red lines. Right: ∇[Fe/H](τ) versus $\mathrm{Range}\stackrel{\sim }{[\mathrm{Fe}/\mathrm{H}]}(\mathrm{age})$ for the same galaxy, with each point representing a different age bin (width 1 Gyr). $\mathrm{Range}\stackrel{\sim }{[\mathrm{Fe}/\mathrm{H}]}(\mathrm{age})$ is determined by looking at the range in [Fe/H] across all age bins, and then normalized so the minimum range is given a value of 0, and the maximum range is given a value of 1. The correlation between $\mathrm{Range}\stackrel{\sim }{[\mathrm{Fe}/\mathrm{H}]}(\mathrm{age})$ and ∇[Fe/H](τ) for this galaxy is −0.87, indicating the evolution of the gradient can be recovered well from $\mathrm{Range}\stackrel{\sim }{[\mathrm{Fe}/\mathrm{H}]}(\mathrm{age})$.

The existence of the correlation between metallicity range and gradient allows for the ability to recover the evolution of the metallicity gradient with cosmic time directly from the data with minimal assumptions. This finding seems natural; the steeper the gradient, the larger the range in [Fe/H]. In our TNG50 sample of Milky Way and Andromeda-like galaxies, we find that most galaxies show a decent correlation between the metallicity gradient and [Fe/H] range (Figure 3), with a median correlation of −0.69. Here, a correlation of −1 implies that ∇[Fe/H](τ) can be perfectly recovered from the scatter in [Fe/H] where the increase in scatter implies a steeper metallicity gradient (see the right panel of Figure 2 as an example of a correlation of −0.86), a correlation of 0 indicates that the scatter in [Fe/H] has no relation with the evolution of the metallicity gradient, and a correlation of 1 implies that the gradient steepens as the range decreases. We also find that this result depends on bar strength, where the stronger barred galaxies have a better ability at recovering ∇[Fe/H](τ) from $\mathrm{Range}\stackrel{\sim }{[\mathrm{Fe}/\mathrm{H}]}(\mathrm{age})$ than weaker barred galaxies.

Fig. 3. Density distribution of the correlation between $\mathrm{Range}\stackrel{\sim }{[\mathrm{Fe}/\mathrm{H}]}(\mathrm{age})$ and ∇[Fe/H](τ) for 179 galaxies (black). These galaxies are chosen such that they have a linear gradient for (at least) their most recent 4 Gyr of star formation. Most galaxies exhibit a reasonable ability to have their ∇[Fe/H](τ) evolution recovered from $\mathrm{Range}\stackrel{\sim }{[\mathrm{Fe}/\mathrm{H}]}(\mathrm{age})$ (correlation ∼ − 1). This ability has a correlation with bar strength (colored lines), where stronger barred galaxies have a better recovery rate than weaker ones.

While most of the sample follows the natural assumption that the metallicity range and gradient are negatively correlated, there are a few galaxies that show little to no correlation between the two variables. There are even some galaxies that show the exact opposite of what is expected; that is, $\mathrm{Range}\stackrel{\sim }{[\mathrm{Fe}/\mathrm{H}]}(\mathrm{age})$ and ∇[Fe/H](τ) are positively correlated, indicating that a steeper gradient has a smaller range in [Fe/H]. These discrepancies are due to an underlying assumption in the ability to connect the scatter in [Fe/H] at a given age to the metallicity gradient at that lookback time – the star-forming region of the disk must be the same radial range across cosmic time. This assumption is not necessarily true in our TNG50 sample (see left panel of Figure 4), though it is crucial because the method claims that a steeper radial metallicity slope creates a larger spread in [Fe/H] than a weaker slope (i.e. a negative correlation). However, if the radial range of star formation (RRSF) was larger at the time of the weaker slope, the scatter in [Fe/H] could be larger. The right panel of Figure 4 illustrates the role that the difference in RRSF plays in our ability to connect $\mathrm{Range}\stackrel{\sim }{[\mathrm{Fe}/\mathrm{H}]}(\mathrm{age})$ and ∇[Fe/H](τ). There is a trend that, on average, the ability to recover the metallicity gradient from the range in [Fe/H] decreases as the RRSF varies drastically over lookback time for a galaxy. There also appears to be a mass dependence, with more massive galaxies showing a better ability to recover ∇[Fe/H](τ) from $\mathrm{Range}\stackrel{\sim }{[\mathrm{Fe}/\mathrm{H}]}(\mathrm{age})$ (similar to what is seen in the NIHAO simulation; Lu et al. 2024).

Fig. 4. Left: Star-forming region for 2 example galaxies (subhalos 506720 and 539333). The shaded region at a specific age indicates the radial range of stars forming during that time. The star-forming region of some galaxies (such as 539333) can vary in width of many kpc, while other galaxies (such as 506720) form stars in a much more consistent range. Right: Correlation between the scatter in [Fe/H] and the metallicity birth gradient as a function of how much the radial range of the star formation (RRSF) varies over time (defined as the difference between the maximum and minimum range stars are formed in the disk over time), colored by the galaxy’s stellar mass within 30 kpc. A value of 0 indicates that the RRSF was constant over time (albeit potentially covering different radii such as subhalo 506720 in the left panel), while a larger value indicates that the galaxy formed stars over a few kpc some Gyr, and over a much larger region other Gyr (such as subhalo 539333 in the left panel). Overall, the more the range of the star-forming region varies over time, the less likely ∇[Fe/H](τ) can be recovered from $\mathrm{Range}\stackrel{\sim }{[\mathrm{Fe}/\mathrm{H}]}(\mathrm{age})$. There is also a gradient with stellar mass, with higher mass galaxies showing a better ability to have their metallicity gradients recovered from $\mathrm{Range}\stackrel{\sim }{[\mathrm{Fe}/\mathrm{H}]}(\mathrm{age})$. The 4 highlighted galaxies in red boxes are examined in detail in Figure 5.

To illustrate how the $\mathrm{Range}\stackrel{\sim }{[\mathrm{Fe}/\mathrm{H}]}(\mathrm{age})$–gradient correlation can be affected by a non-uniform RRSF over time, Figure 5 shows the relative strength of the [Fe/H] gradient versus the RRSF, with the correlation between $\mathrm{Range}\stackrel{\sim }{[\mathrm{Fe}/\mathrm{H}]}(\mathrm{age})$ and the metallicity gradient shown in the bottom left corner. The four galaxies chosen are highlighted in red in the right panel of Figure 4, and show examples of when the method works regardless of variation in RRSF over time (top row of Figure 5) and when the method does not work (bottom row of Figure 5). Specifically, the top left panel depicts a galaxy whose RRSF stays relatively constant with time. Thus, as the [Fe/H] gradient evolves, there is a direct correlation with the range in [Fe/H]. The top right panel also shows a galaxy with a strong ability to recover the [Fe/H] gradient from the range, but the RRSF increases by ∼8 kpc throughout its lifetime. Here, the metallicity gradient is weakest when the RRSF is smallest, and the gradient steepens with increasing RRSF, further enforcing a strong correlation between $\mathrm{Range}\stackrel{\sim }{[\mathrm{Fe}/\mathrm{H}]}(\mathrm{age})$ and ∇[Fe/H](τ).

Fig. 5. Radial range of star-forming region (RRSF) versus the relative strength of the metallicity gradient for different ages, for the four highlighted galaxies in Figure 4. The correlation between ∇[Fe/H](τ) and $\mathrm{Range}\stackrel{\sim }{[\mathrm{Fe}/\mathrm{H}]}(\mathrm{age})$ is given in the bottom left of each panel, below the subhalo ID. The top panels depict situations in which the evolution of the metallicity gradient can be recovered from the scatter in [Fe/H], while the bottom row illustrates two examples of galaxies where the gradient cannot be recovered in this fashion due to the gradient being steepest during a time of a lower star-forming region.

There are some outliers in Figure 4 that show ∇[Fe/H](τ) cannot be recovered using the method of Lu et al. (2022a) even with a minimally varying RRSF. The bottom row of Figure 5 depicts two examples where the method does not work; one of a galaxy with essentially no relation between ∇[Fe/H](τ) and $\mathrm{Range}\stackrel{\sim }{[\mathrm{Fe}/\mathrm{H}]}(\mathrm{age})$ (bottom left) and one with a positive correlation between the two variables (bottom right). For the bottom left subhalo, even though the RRSF only varies by a small amount (∼2 kpc) across cosmic time, the gradient evolves such that the steepest gradient happens when the new stars are forming in a smaller R_birth range than when the gradient is weakest. A more extreme example is shown in the bottom right panel, where, in this case, the RRSF varies a lot. Overall, for TNG50 Milky Way and Andromeda-like galaxies whose RRSF varies by < 8 kpc throughout the disk’s lifetime, we find an excellent ability to recover the time evolution of the metallicity gradient from the range in [Fe/H] over age (median correlation of −0.9). When the RRSF varies more throughout the galaxy’s lifetime (> 8 kpc), the median correlation drops down to −0.51.

4.3. How galaxies form affect ability to recover gradient

Section 4.2 illustrates that the evolution of the [Fe/H] gradient can be recovered from the range in [Fe/H] across age for stronger-barred galaxies, while weaker-barred galaxies exhibit little to no correlation. The lack of correlation between the [Fe/H] gradient and range is driven by the width of the star forming region varying across time (Fig. 4). The top panel of Figure 6 shows how the width of the star forming region changes across time for different barred galaxies. Overall, the weaker-barred galaxies see an increase in the width of the star forming region, where these galaxies begin forming stars on the inner few kpc, and then grow over time until they are forming stars throughout essentially the full disk. This is also illustrated in the middle and bottom panels of the same figure, with the middle panel depicting the growth of the disk while the bottom panel shows that weaker-barred galaxies consistently form stars in the inner regions.

Fig. 6. Top: Growth of the regions where stars form in the galactic disk over time, separated by bar strength at redshift 0. Weaker-barred galaxies start with stars forming on less than 3 kpc in the center of the disk, and then grow until stars are forming almost throughout the entirety of the disk. Stronger-barred galaxies show only minimal growth in their radial range of star-forming region (RRSF) Middle: Upper radial limit of where stars are forming throughout the galaxy at a given moment. Each group of galaxies shows the expected inside-out growth. Bottom: Lower radial limit of where stars are forming in the galaxy over time. Weaker-barred galaxies tend to form stars consistently in their inner regions, whereas stronger-barred galaxies stop forming stars in the inner few kpc. In each panel, the lines correspond to the median of each group and the shaded region represents the 25−75% quantiles. To best illustrate the growth of the disk over time, we use age bins of 0.2 Gyr.

Stronger-barred galaxies, on the other hand, show a different evolution. While they appear to have a similar inside-out growth as other galaxies (middle panel Figure 6), stronger-barred galaxies exhibit quenching in the bar region starting about 6−7 Gyr ago. Over time, the bar and lower radius of star formation increase (bottom panel Figure 6), causing the RRSF to remain relatively constant during the disk’s evolution, and allowing for a better recovery of the metallicity gradient from the range in [Fe/H] (right panel of Figure 4). This halt in star formation in the inner regions of barred galaxies is seen in other simulations (e.g. Khoperskov et al. 2018) and observational data (Géron et al. 2024).

5. Results II: recovering the evolution of [Fe/H] at the galactic center

Now that we have established when the evolution of the metallicity gradient can be recovered from the range in [Fe/H] across age, we move to the other key component in estimating the metallicity profile over time and radius – the evolution of [Fe/H] in the inner few kpc. In this section, we investigated the ability to recover the evolution of central metallicity (Section 5.1) and directly tested the assumption that [Fe/H](0, τ) is monotonically increasing over time (Section 5.2). This section also argues for the use of a projected central metallicity when recovering the galactic disk metallicity profile.

5.1. Ability to recover the central metallicity with time

As discussed in Section 3, estimating birth radii requires knowledge of [Fe/H](R_birth = 0 kpc, τ). In recent methods, this quantity is estimated from the upper envelope of the age-metallicity relation for older stars, and extrapolated for younger (≲8 Gyr) stars. Since current samples do not have young stars that were formed in the Galactic center, the boundary condition for this extrapolation is that [Fe/H](R_birth = 0 kpc, τ = 0 Gyr) is estimated from the metallicity in the solar neighborhood and the present day gradient of the youngest stars. Figure 1 illustrated that on average, the galaxies of TNG50 have a fairly linear gradient. However, as these gradients are not completely linear, we wish to test our ability in recovering [Fe/H](R_birth = 0 kpc, τ). While the method of Lu et al. (2022a) only uses this technique to recover the boundary condition for the present day, here we investigate it for all age bins that the galaxy has an on average linear gradient to gather a more complete picture.

Figure 7 illustrates an example of how [Fe/H](R_birth = 0 kpc, τ) is estimated for a given period of time. Using the metallicity of a given radius (in this case 8.5 kpc; black square) and the metallicity gradient at that time, we can estimate the projected [Fe/H](R_birth = 0 kpc) ($\widehat{[\mathrm{Fe}/\mathrm{H}]}(R_{\mathrm{birth}}$ = 0 kpc); black star). The true central metallicity is taken as the median [Fe/H](R_birth < 1 kpc) of disk stars, and is shown as the red circle. In this example, the estimated central metallicity underpredicts the truth by nearly 0.4 dex due to the large spike in [Fe/H] in the galactic center. Figure 8 shows the density distribution of the difference between the true and projected central metallicities across all times and galaxies where the galaxies have stars forming at R_birth < 1 kpc, with the top panel separated by bins of lookback time and the bottom panel separated by bar strength at redshift 0. The differences between the true and projected central metallicity can be quite large (up to 0.5 dex), with a bias towards negative differences. This bias indicates that the projected central metallicity typically underestimates the true central metallicity, which is precisely what was witnessed in Figure 7. The offset is more noticeable for larger lookback times and weaker barred galaxies, though the underestimation exists for nearly all instances.

Fig. 7. Illustrating the recovery of [Fe/H](Rbirth = 0 kpc) for subhalo 553837 for stellar particles 7−8 Gyr ago. The projected central metallicity ($\widehat{[\mathrm{Fe}/\mathrm{H}]}(R_{\mathrm{birth}}$ = 0 kpc); black star) was estimated from the metallicity at a given radius (8.5 kpc; black square) and the metallicity gradient at the current time. The red circle represents the true central metallicity [Fe/H](Rbirth = 0 kpc), taken as the median [Fe/H](Rbirth < 1 kpc). For this scenario, the central metallicity cannot be estimated accurately assuming a linear gradient due to the spike in [Fe/H] in the inner kpc. This galaxy has a weak bar (bar strength of 0.15) and a mass of 3.5 × 1010 M⊙.

Fig. 8. Density distribution of the difference between inferred and true [Fe/H] at galactic center separated by top: lookback time τ and bottom: bar strength. The projected [Fe/H](R = 0 kpc, τ) is estimated from the [Fe/H] value at a given radius and the current [Fe/H] gradient for a given age of width 1 Gyr (see Figure 7). As a whole, the value of [Fe/H](R = 0 kpc, τ) is underestimated by ∼0.2 dex, with larger differences for weaker barred galaxies and larger lookback times. The bias towards negative values (i.e. [Fe/H](R = 0 kpc, τ) is underestimated) is predominantly due to the non-linear spike in [Fe/H] in the inner kpc (see Figure 7 for an example).

The consistency in the projected central metallicity underestimating the true central metallicity illustrates that the spike in the galactic central region’s [Fe/H] is a stable feature, at least in TNG50 Milky Way and Andromeda analogues. While the difference can be up to ∼ − 0.5 dex, this underestimation will only affect the recovery of the metallicity evolution for the inner kpc. The remainder of the disk’s [Fe/H] profile is fairly well estimated (such as in Figure 7). In terms of estimating birth radii, the stars in this [Fe/H]-spike region will be given negative birth radii, and thus their birth radii are not able to be recovered with current assumptions. Although the centres of the TNG50 galaxies may not be perfectly resolved, keeping in mind the mass resolution of 8 × 10⁴ M_⊙, the simulations are still likely to capture differences in physical conditions for star formation and its efficiency relative to the rest of the disk (Querejeta et al. 2021; Pessa et al. 2023). This may result in starbursts, especially at higher redshifts when mergers and interactions are more abundant (Renaud et al. 2019). For instance, the higher pressure in the galactic nuclei is likely to result in a more efficient star formation where the released enriched gas, while being captured in the potential well, recirculated on a shorter time scale, thus resulting in a rapid but local enhancement in the metallicity of newly formed stars. However, such a process can be counteracted by the AGN feedback or other types of star formation quenching; though, it is not the case in many galaxies since we observe continuous star formation in their centres.

5.2. Projected [Fe/H] evolution at the galactic center

A main assumption in the birth radii method used in Minchev et al. (2018), Lu et al. (2022a), Ratcliffe et al. (2023, 2024) is that the metallicity at the Galactic center ([Fe/H](R_birth = 0 kpc, τ)) is monotonically increasing. However, this assumption can be violated in simulations, where gas-rich mergers cause a dilution in the metallicity of stars forming throughout the host galaxy (e.g. Buck 2020; Agertz et al. 2021; Annem & Khoperskov 2024; Renaud et al. 2024). This assumption also proves false in the majority of the galaxies in our TNG50 Milky Way and Andromeda-like sample, where we find over 90% of the galaxies in our sample show a non-monotonic evolution in [Fe/H] in the galactic center (R_birth < 1 kpc) to redshift 0. The left schematic of Figure 9 illustrates the evolution of the true central metallicity for each galaxy. Starting with 183 galaxies (15 galaxies have too few ages with R_birth < 1 kpc), we find 15 do not show a depletion in [Fe/H](R_birth = 0 kpc, τ) down to redshift 0, while 45 galaxies show at least one major depletion in [Fe/H] in their galactic centers. The remaining 123 galaxies exhibit quenching at some point in their central kpc, and no stars are formed with R_birth < 1 kpc after. The mechanisms of the star formation quenching in the galactic center can be diverse, mainly including morphological (Martig et al. 2009; Renaud et al. 2015; Khoperskov et al. 2018) and AGN-feedback (Nelson et al. 2018, 2019b; Weinberger et al. 2018), whose relative importance across the sample of our galaxies requires further investigation and is outside the scope of this work.

Fig. 9. Schematic illustrating the distribution of galaxies and how their left: [Fe/H](R = 0 kpc, τ) and right: $\widehat{[\mathrm{Fe}/\mathrm{H}]}(R_{\mathrm{birth}}$ = 0 kpc) evolves over time. The evolution of the true central metallicity (left) was estimated by tracing how [Fe/H] in the central kpc of the galactic disk sample enriches over time. Of the 183 galaxies that have enough disk stars born within the inner kpc, only 15 of them are found to have their centers predominantly continuously enrich in [Fe/H] over cosmic time, while 45 of the galaxies are found to have at least one major depletion event in their galactic centers. The remaining 123 galaxies had star formation quench before redshift 0 in their center. An example of each of these scenarios is given, where the density of the age–metallicity relation for our defined disk sample (grey) and stellar particles with Rbirth < 1 kpc (blue) is shown. The evolution of the projected central metallicity (right) was determined by tracing the projected [Fe/H] at Rbirth = 0 kpc using the metallicity at a given radius and the current metallicity gradient (see Figure 7 for an example). $\widehat{[\mathrm{Fe}/\mathrm{H}]}(R_{\mathrm{birth}}$ = 0 kpc) traces the central metallicity assuming a linear metallicity gradient, and acts as an anchor point for the disk’s gradient, regardless if the true central metallicity spikes in the galactic center. Of the 178 galaxies with enough disk stars to appropriately trace the time evolution of $\widehat{[\mathrm{Fe}/\mathrm{H}]}(R_{\mathrm{birth}}$ = 0 kpc), we find 81 galaxies show a continuous enrichment at Rbirth = 0 kpc, while 97 have at least one major depletion event. Examples of these two scenarios are given next to the respective outcome, with the blue line indicating the projected central metallicity, estimated using age bins of 0.5 Gyr. The central metallicity of the Milky Way is assumed to continuously increase in current methods, with potentially only minor deviations from this assumption (see Section 6.1.1).

Figure 9 illustrates that most galaxies have star formation halted in their galactic centers sometime during their evolution. This shows that tracing the upper boundary of [Fe/H] in the age–metallicity plane across all ages will rarely provide the evolution of [Fe/H](R_birth = 0 kpc, τ), even with the most complete sample. Additionally, the previous section (Section 5.1) discussed how the metallicity at the galactic center is not necessarily representative of the metallicity profile in the remainder of the disk, at least for our TNG50 Milky Way and Andromeda-like sample. Thus, we propose not to consider the true central metallicity when deriving the disk’s metallicity evolution with time, and hence birth radii, but rather use the projected metallicity $\widehat{[\mathrm{Fe}/\mathrm{H}]}(R_{\mathrm{birth}}$ = 0 kpc, τ). Focusing on the evolution of the projected metallicity (right schematic in Figure 9), 46% of the galaxies show a continuously enriching projected central metallicity, and 54% of the galaxies exhibit at least one major dilution.

Figure 9 informs us that many galaxies experience a dilution in their inner kpc and projected $\widehat{[\mathrm{Fe}/\mathrm{H}]}(R_{\mathrm{birth}}$ = 0 kpc, τ), and break the assumption of a monotonically enriching galactic center. To investigate potential reasons why, we show the distribution of times between the depletion in $\widehat{[\mathrm{Fe}/\mathrm{H}]}(R_{\mathrm{birth}}$ = 0 kpc, τ) and when the galaxy experienced its last major (blue) and minor/major (orange) in Figure 10. In the sample of 198 galaxies, we find 131 instances of the central metallicity decreasing by at least 0.05 dex in 0.5 Gyr. This indicates that some galaxies experienced more than one depletion in their central [Fe/H]. Of these drops, 79% of them had a merger with a gas mass ratio > 10% in the 2 Gyr prior to the depletion, and 58% of these instances had a massive merger (gas mass ratio > 25%) within the previous 2 Gyr. This shows that there is a correlation between mergers and depletion in $\widehat{[\mathrm{Fe}/\mathrm{H}]}(R_{\mathrm{birth}}$ = 0 kpc, τ).

Fig. 10. Distribution of times between the depletion event in $\widehat{[\mathrm{Fe}/\mathrm{H}]}(R_{\mathrm{birth}}$ = 0 kpc, τ) and the last major (blue; gas mass ratio > 25%) and minor/major (orange; gas mass ratio > 10%) mergers. A depletion event is defined as the projected central metallicity decreasing by at least 0.05 dex in 0.5 Gyr, smoothed over 2.5 Gyr. The majority of the time that the projected metallicity in the galactic central region decreases, the galaxy experienced a merger of at least 10% gas mass fraction within the previous 2 Gyr. This suggests that larger mergers potentially affect the evolution of the metallicity of the galactic center.

Figure 10 answers the question of if there was a merger prior to a depletion event. However, this does not necessarily mean that mergers always cause dilutions in $\widehat{[\mathrm{Fe}/\mathrm{H}]}(R_{\mathrm{birth}}$ = 0 kpc, τ). To answer what happens to the projected central metallicity after a merger, Figure 11 shows the change in $\widehat{[\mathrm{Fe}/\mathrm{H}]}(R_{\mathrm{birth}}$ = 0 kpc, τ) as a function of time after a merger event for each galaxy after a merger, until either another merger happens or redshift 0. While some galactic centers dilute in projected [Fe/H] immediately after a merger, there is the general trend that it can take up to a few Gyr before $\widehat{[\mathrm{Fe}/\mathrm{H}]}(R_{\mathrm{birth}}$ = 0 kpc, τ) dilutes. This result suggests that mergers do not necessarily cause the projected central metallicity to dilute, but if there is a dilution, it is not instantaneous. It also shows that mergers can slow down enrichment in the projected central metallicity, even if there is not a dilution.

Fig. 11. Density of the change in $\widehat{[\mathrm{Fe}/\mathrm{H}]}(R_{\mathrm{birth}}$ = 0 kpc, τ) as a function of time after a (left) massive or (right) minor merger for all galaxies after all massive/minor mergers, until either another massive/minor merger occurs or no more stars form in the central 1 kpc. A positive change in [Fe/H] indicates that the projected metallicity in the central 1 kpc increased, while a negative value indicates that the projected [Fe/H] in the central kpc decreased. The projected metallicity at the galactic center does not necessarily decrease immediately after a merger, and a dilution in $\widehat{[\mathrm{Fe}/\mathrm{H}]}(R_{\mathrm{birth}}$ = 0 kpc, τ) can take up to a few Gyr before diluting the projected galactic center in [Fe/H]. Even if $\widehat{[\mathrm{Fe}/\mathrm{H}]}(R_{\mathrm{birth}}$ = 0 kpc, τ) does not dilute, there is a general trend that [Fe/H] enrichment slows down.

In addition to showing that mergers may not cause an instantaneous decrease in the projected central metallicity, Figure 11 also shows that $\widehat{[\mathrm{Fe}/\mathrm{H}]}(R_{\mathrm{birth}}$ = 0 kpc, τ) can continue to increase after a merger. In our investigations, we find no overarching relationship between merger properties and their effect on the projected metallicity. For instance, when the merger brings in pre-enriched gas, the overall metallicity of the disk can increase, causing an increase in the projected metallicity (top panels of Figure 12). However, a merger with pre-enriched gas can also cause the projected central metallicity to decrease (bottom panels of Figure 12). This is due to the enriched gas causing the [Fe/H] gradient to weaken, and therefore the extension of the disk to the center is thus lowered. If the merger contains pristine gas, the overall metallicity of the disk can decrease (as seen in e.g. VINTERGATAN; Renaud et al. 2024) and/or the gradient of the disk can steepen (as discussed in Buck et al. 2023) and therefore create an increase in the projected central metallicity. A detailed paper providing an in depth investigation into what types of mergers cause a depletion or continual increase would provide better insights into the conclusions we draw here.

Fig. 12. Two examples of the evolution of $\widehat{[\mathrm{Fe}/\mathrm{H}]}(R_{\mathrm{birth}}$ = 0 kpc, τ) after a merger with pre-enriched gas. Due to the metal-rich spike in the galactic center (Figure 7), we remove the inner kpc for visualization purposes. The subhalo in the top row (subhalo 554523) had a major merger at lookback time ∼9 Gyr ago which brought in more metal-rich gas and increased the metallicity of the entire disk. This allowed for the projected central metallicity (marked with the red “x”) to also increase. The subhalo in the bottom row (subhalo 514829) experienced a metal-rich merger ∼8 Gyr ago. However, instead of increasing the metallicity throughout the entire disk like in the top row, here only the outer disk had an increase in [Fe/H]. This thus weakened the metallicity gradient, and in turn lowered the projected central metallicity. These examples show that even something as trivial as the metallicity of the merger compared to the main galaxy cannot begin to disentangle when $\widehat{[\mathrm{Fe}/\mathrm{H}]}(R_{\mathrm{birth}}$ = 0 kpc, τ) increases or decreases after a merger.

6. Discussion

6.1. Probability of success in recovering birth radii for disk stars in the Milky Way

While hydrodynamical cosmological simulations and Galactic chemical evolution models have made great advances in understanding the Milky Way’s formation and evolution, there are still discrepancies that only exploring the Milky Way itself will answer. Knowledge of stellar birth radii in the Milky Way disk provides a useful avenue in inferring the disk’s evolution, from understanding the strength of migration (e.g. Frankel et al. 2018) and mergers on the Galactic disk (Lu et al. 2022a), inferring the chemical evolution of the disk (Ratcliffe et al. 2023, 2024; Wang et al. 2024), and constraining chemo-kinematic relations (Minchev et al. 2018). However, we cannot directly measure R_birth, and need to infer it using some method (Minchev et al. 2018; Frankel et al. 2018; Lu et al. 2022a; Wang et al. 2024).

The focus of this paper is to investigate the recent method from Lu et al. (2022a) (whose assumptions are quite similar to previous methods and necessary for R_birth recovery) with TNG50 Milky Way and Andromeda-like galaxies, and understand its success in the context of the Milky Way itself. The method uses Equation (2), where ∇[Fe/H](τ) is estimated from the range in [Fe/H] over different age bins (Section 4.2) and [Fe/H](R = 0 kpc, τ) is estimated from the upper envelope of older stars (when the sample is expected to contain stars that formed in the inner few kpc of the Galaxy) and modeled as continuously enriching in [Fe/H] until the present day, with the boundary condition estimated using the current metallicity and gradient of the youngest stars in the solar neighborhood (Section 5.1). Overall, we find that the individual aspects of this method seem like they can be reasonably extended to the Milky Way disk. Most galaxies in the sample show that the metallicity gradient can be recovered from the range in [Fe/H] across age (Figure 3). The galaxies that fail this condition tend to have weaker bars compared to the stronger one of the Milky Way (at least in the context of TNG50; Khoperskov et al. 2023b) and have large variations in the range of their star forming region over time (> 8 kpc). While the median correlation between $\mathrm{Range}\stackrel{\sim }{[\mathrm{Fe}/\mathrm{H}]}(\mathrm{age})$ and ∇[Fe/H](τ) is −0.85 for galaxies whose star-forming region can vary in size by up to 8 kpc, a more accurate method would account for the growth of the Galactic disk (e.g. Mollá et al. 2019).

While the metallicity gradient provides the relative [Fe/H] across the disk, [Fe/H](R = 0 kpc, τ) acts as an anchor in determining the precise [Fe/H] as a function of radius, even if there is no star formation in the galactic center over the previous few Gyr. Section 5 discussed the importance of using the projected central metallicity ($\widehat{[\mathrm{Fe}/\mathrm{H}]}(R_{\mathrm{birth}}$ = 0 kpc, τ)) due to the possibly non-representative [Fe/H] profile in the inner few kpc. Section 5.2 suggests that the assumption that $\widehat{[\mathrm{Fe}/\mathrm{H}]}(R_{\mathrm{birth}}$ = 0 kpc, τ) is continuously enriching in metallicity is potentially violated for the average Milky Way and Andromeda-like galaxy, since we find just over 50% of the galaxies in our sample have at least one dilution in the projected [Fe/H] in their galactic centers at some point during their evolution after the disk formed. However, the vast majority of the times when the galactic centers see a depletion in projected [Fe/H] are within 2 Gyr of a merger with gas mass ratio > 10% (Figure 10). This suggests that the Milky Way’s Galactic center probably did not see a decrease in [Fe/H] with time due to Sagittarius. Though, it does raise the question of the validity of $\widehat{[\mathrm{Fe}/\mathrm{H}]}(R_{\mathrm{birth}}$ = 0 kpc, τ) monotonically increasing after the GSE merger event.

6.1.1. Incorrect projected central metallicity assumption gives non-physical birth radii

An additional constraint in the R_birth methods from Minchev et al. (2018), Lu et al. (2022a) is that the birth radii distributions for mono-age populations in the solar neighborhood need to be physically meaningful (see Figure A.3 in Lu et al. 2022a). In particular, the youngest stars in the solar neighborhood are expected to be locally born, while the oldest stars are expected to be born more inwards (Minchev et al. 2013; Netopil et al. 2022; Agertz et al. 2021; Carrillo et al. 2023). Therefore, a natural test to determine the effect of incorrectly assuming a monotonically increasing $\widehat{[\mathrm{Fe}/\mathrm{H}]}(R_{\mathrm{birth}}$ = 0 kpc, τ) is to see how the R_birth distributions of mono-age populations of stellar particles near 8 kpc vary with this assumption. Figure 13 shows this test for one of our TNG50 Milky Way and Andromeda-like galaxies (subhlao 560751). This galaxy is like the Milky Way in many ways: it has a stronger bar strength (0.61), M_*(< 30 kpc) = 4.12 × 10¹⁰ M_⊙, and it appears to have no stars forming in the inner few kpc for the last ∼5 Gyr. The projected metallicity in the galactic center of this galaxy does not monotonically increase, as illustrated by the black line in the age–metallicity relation shown in the left panel of Figure 13. The red line in the same panel shows what central metallicity evolution is assumed for the Milky Way in Lu et al. (2022a), Ratcliffe et al. (2023, 2024); the upper [Fe/H] envelope for older ages and a log function for younger ages with the boundary condition that the $\widehat{[\mathrm{Fe}/\mathrm{H}]}(R_{\mathrm{birth}}$ = 0 kpc, τ) is estimated from the present day gradient of the youngest stars and the metallicity at a given radius. We estimated R_birth in two different ways for this galaxy – using the true, declining $\widehat{[\mathrm{Fe}/\mathrm{H}]}(R_{\mathrm{birth}}$ = 0 kpc, τ) (black line from left panel of Figure 13) which we call Inferred_decR_birth and using the estimated, inclining $\widehat{[\mathrm{Fe}/\mathrm{H}]}(R_{\mathrm{birth}}$ = 0 kpc, τ) (red line from same figure) which we call Inferred_incR_birth. For both ways, we used the gradient estimated from $\mathrm{Range}\stackrel{\sim }{[\mathrm{Fe}/\mathrm{H}]}(\mathrm{age})$ in Equation (2). To make the best connections between the simulation and the Milky Way, we estimated R_birth after redrawing new age and [Fe/H] parameters for each stellar particle by using uncertainties of 15% and 0.05 dex respectively. These errors are on par with the measurement uncertainties of large spectroscopic surveys (Jönsson et al. 2020; Xiang et al. 2019; Queiroz et al. 2023; Anders et al. 2023).

Fig. 13. Ability to recover Rbirth for subhalo 560751 with age and [Fe/H] uncertainties added (15% and 0.05 dex respectively). Left: The age–metallicity relation for this subhalo with the true evolution of $\widehat{[\mathrm{Fe}/\mathrm{H}]}(R_{\mathrm{birth}}$ = 0 kpc, τ) (black line) and estimated evolution assuming the projected galactic center is always increasing in [Fe/H] (red line). The second, third, and right panels show the Rbirth distributions for mono-age populations of solar neighborhood (7.75 ≤ R ≤ 8.25) stellar particles for the true Rbirth, Rbirth inferred using the true $\widehat{[\mathrm{Fe}/\mathrm{H}]}(R_{\mathrm{birth}}$ = 0 kpc, τ) (InferreddecRbirth), and Rbirth inferred using the monotonically increasing $\widehat{[\mathrm{Fe}/\mathrm{H}]}(R_{\mathrm{birth}}$ = 0 kpc, τ) (InferredincRbirth) respectively. Inappropriately assuming a monotonically increasing metallicity abundance in the galactic center provides Rbirth distributions that are non-physical and do not match the true distributions, whereas using the true $\widehat{[\mathrm{Fe}/\mathrm{H}]}(R_{\mathrm{birth}}$ = 0 kpc, τ) evolution provides a more accurate trend in how Rbirth decreases with increasing age. Examining the Rbirth distributions of mono-age populations in the solar neighborhood is thus a necessity in constraining $\widehat{[\mathrm{Fe}/\mathrm{H}]}(R_{\mathrm{birth}}$ = 0 kpc, τ) evolution.

The other panels of Figure 13 illustrate the R_birth distributions of mono-age populations for stellar particles within 7.75 ≤ R ≤ 8.25 kpc for the true R_birth (second panel), Inferred_decR_birth (third panel), and Inferred_incR_birth (right panel). The distributions of the true R_birth show clear signs of radial migration and inside-out formation with the oldest stars forming R_birth < 5 kpc and the younger stars forming in the 8 kpc region. Recovering R_birth using the estimated gradient and declining $\widehat{[\mathrm{Fe}/\mathrm{H}]}(R_{\mathrm{birth}}$ = 0 kpc, τ) (Inferred_decR_birth) shows similar trends with age. Assuming a monotonically increasing $\widehat{[\mathrm{Fe}/\mathrm{H}]}(R_{\mathrm{birth}}$ = 0 kpc, τ) provides incorrect distributions, where now the stars forming 4 − 8 Gyr ago are estimated to have formed beyond 8 kpc, with a large percentage supposedly forming beyond 15 kpc. Additionally, the youngest (≤2 Gyr) are found to be born at ∼10 kpc, implying these stars migrated inwards 2 kpc in less than 2 Gyr. Similarly wrong results still exist even when we change the strength of the gradient (see Figure B.5 in the appendix), which has been shown to affect the R_birth distributions of mono-age populations (Lu et al. 2022a). This finding holds good news for extending the R_birth method to the Milky Way disk. The R_birth distributions of mono-age populations for solar neighborhood stars in the Milky Way have been found to show a clear decrease in R_birth with increasing age (see Fig. 6 in Lu et al. 2022a and Fig. 3 in Ratcliffe et al. 2023), much like the distributions found in the TNG50 galaxy for the R_birth estimated using the decreasing $\widehat{[\mathrm{Fe}/\mathrm{H}]}(R_{\mathrm{birth}}$ = 0 kpc, τ); the distributions for Milky Way disk stars does not resemble the distributions found when incorrectly assuming the galactic center’s metallicity monotonically increases. This finding strongly suggests that the metallicity evolution of the Milky Way’s inner few kpc continuously enriches in [Fe/H], with potentially only minor deviations.

6.1.2. Assuming Milky Way-like scatter about the metallicity gradient

The third panel of Figure 13 shows that even though we are able to correctly recover the trend between age and R_birth for solar neighborhood stars, the recovered distributions are wider than the truth. This is in part due to the scatter seen in the simulation about the linear metallicity gradient, which we find to be 0.14 dex across all ages for subhalo 560751. This scatter can be explained by strong radial migration, which has been shown to be unrealistically strong in TNG50 simulated galaxies (Wang & Lilly 2023), and may be due to the rapid migration from feedback in simulations (El-Badry et al. 2016). To account for this effect, we now adjust the dispersion about the linear birth metallicity gradient for every age bin (1 Gyr) in TNG50 galaxy 560751 to be on the same order as that seen in H II regions of the Milky Way (0.07 dex; Esteban et al. 2017; Arellano-Córdova et al. 2021). Figure 14 shows how much more successful the R_birth method is when the dispersion is taken to be more Milky Way-like. The left panel shows the inferred vs true birth radii when R_birth is estimated using Equation (2) with the estimated ∇[Fe/H](τ) and the decreasing $\widehat{[\mathrm{Fe}/\mathrm{H}]}(R_{\mathrm{birth}}$ = 0 kpc, τ) (black line in Figure 13). Even in the regime of extra migration from feedback, we are able to recover R_birth to within 1.34 kpc. When the [Fe/H] dispersion is adjusted to follow that of the Milky Way, the ability to recover R_birth increases to within 0.67 kpc. This value is on-par with results using higher resolution simulations in recovering R_birth (Lu et al. 2022a). The positive bias implies that the method overestimates the birth radius of a star on average by 0.42 kpc, which we find shrinks with increasing precision in age and [Fe/H].

Fig. 14. Inferred Rbirth versus true Rbirth for subhalo 560751 using the method from Lu et al. (2022a) with age and [Fe/H] uncertainties added of 15% and 0.05 dex respectively. Rbirth are estimated using the true evolution of [Fe/H](R = 0 kpc, τ) (black line in Fig. 13) and ∇[Fe/H](τ) recovered from $\mathrm{Range}\stackrel{\sim }{[\mathrm{Fe}/\mathrm{H}]}(\mathrm{age})$. The difference between the left and right panels is that stellar [Fe/H] values in the TNG50 galaxy were adapted to better match the dispersion about the metallicity gradient as seen in the Milky Way’s H II regions in the right panel. This illustrates the high probability of success in estimating Rbirth in the Milky Way disk.

6.2. Do wiggles happen at same age

When calculating stellar birth radii for APOGEE DR17 (Abdurro’uf et al. 2022; Majewski et al. 2017) red giants, Ratcliffe et al. (2023) found that a given birth radius did not necessarily continuously enrich in [Fe/H] over time (see the right panel of their Fig. 2). These “wiggles” were seen to happen at the same time for each birth radius (in their case, ∼4 and 6 Gyr) with the outer birth radii having the larger fluctuation in [Fe/H]. This structure is caused when [Fe/H](R = 0 kpc, τ) minimally varies when the metallicity gradient fluctuates (i.e. flattens, steepens, and then flattens again). While Section 5.2 showed that most galaxies in our sample do not monotonically increase in [Fe/H] at their center, we still wish to understand if these “wiggles” do in fact happen at the same time for each birth location, or if the outer (or inner) radii feel this effect first.

The left column of Figure 15 shows the age–metallicity relation for mono-R_birth populations for a galaxy from our sample whose mono-R_birth populations do not continuously enrich in [Fe/H] and “wiggle”. Due to the variability in [Fe/H](R = 0 kpc, τ), examining the extent of these wiggles is nontrivial. Therefore, to better match the conditions assumed for the Milky Way disk during the time of gradient fluctuations in Ratcliffe et al. (2023), we also show the metallicity of a mono-R_birth population relative to the metallicity of the inner 2 kpc in the right column of Figure 15. Once [Fe/H](R, τ) is put relative to the metallicity of the inner few kpc, the wiggles become much more apparent. It is easier to see that the wiggles happen at the same time, similar to what what shown in Ratcliffe et al. (2023). Most other galaxies that exhibit wiggles follow similar trends. The wiggle at 4 Gyr may be due to two metal-poor mergers that merged ∼5 − 6 Gyr ago (gas merger ratios of 4% and 10%), while the dramatic increase in [Fe/H] in the outer disk at ∼3 Gyr may be due to the few small metal-rich mergers at ∼4 − 5 Gyr. Similar trends can also be seen in the NIHAO-UHD simulation suite (Buck et al. 2023).

Fig. 15. Left column: Evolution of mono-Rbirth populations in the age–metallicity relation for a TNG50 Milky Way and Andromeda-like galaxy (subhalo 402555) whose mono-Rbirth populations exhibit “wiggling” behavior, where each mono-Rbirth population did not experience a monotonic enrichment in [Fe/H]. Right column: [Fe/H] abundance of mono-Rbirth populations relative to [Fe/H](R < 2 kpc) across age. Looking at the [Fe/H] abundance relative to the inner few kpc at a given age allows the structure of the wiggles to become more clear. Most of the galaxies that exhibit “wiggling” structure in this plane follow similar trends to this galaxy – each Rbirth sees a depletion in [Fe/H] at the same time as also seen using Rbirth methods (see Fig. 2 in Ratcliffe et al. 2023).

7. Conclusions

This work explores in detail the validity of assumptions needed to estimate [Fe/H](R, τ) and R_birth using the sample of TNG50 Milky Way and Andromeda analogues, and discusses the implications of these assumptions in the context of the Milky Way. Knowledge of the evolution of the metallicity profile with radius, in addition to stellar R_birth, opens up a new and powerful avenue to explore the Milky Way’s formation path, and can provide necessary constraints on chemo-dynamical evolution. Our key findings are:

• Most Milky Way and Andromeda-like galaxies in TNG50 have a good correlation between [Fe/H] and R_birth across time (Figure 1). We find that most of the galaxies that do not exhibit a tight correlation either have too shallow of a gradient (which causes a low correlation) or seemingly had their pre-existing disk destroyed by late massive mergers.

• We can recover the time evolution of the metallicity gradient for most TNG50 Milky Way and Andromeda-like galaxies in our sample from the range in [Fe/H] across age, with a dependency on the bar strength at redshift zero (Figure 3). The galaxies that show no correlation between ∇[Fe/H](τ) and $\mathrm{Range}\stackrel{\sim }{[\mathrm{Fe}/\mathrm{H}]}(\mathrm{age})$ have a large variability in the radial range of their star-forming region over time. The galaxies in our sample with weaker bars (bar strength < 0.2) predominantly have larger variations in their star-forming region, while the stronger barred galaxies have less variation on average.

• [Fe/H](R = 0 kpc, τ) is not necessarily representative of the metallicity profile of the remainder of the disk due to, for example, spikes in [Fe/H] in the inner few kpc (Figures 7 and 8). The differences are smaller for smaller lookback times and galaxies with stronger bars. We propose instead that recovering the evolution of the metallicity profile of the galactic disk should use a projected central metallicity ($\widehat{[\mathrm{Fe}/\mathrm{H}]}(R_{\mathrm{birth}}=0$ kpc, τ)) that captures the trend of the disk and ignores non-linearities in the center.

• The evolution of $\widehat{[\mathrm{Fe}/\mathrm{H}]}(R_{\mathrm{birth}}$ = 0 kpc, τ) is non-monotonic for ∼55% of the galaxies in our sample (Figure 9). 79% of these dilutions in [Fe/H] in the galactic center happened within 2 Gyr of a merger with gas mass ratio > 10% (Figure 10).

• Requiring R_birth distributions of mono-age populations to remain physically meaningful constrains the assumption of a monotonically increasing $\widehat{[\mathrm{Fe}/\mathrm{H}]}(R_{\mathrm{birth}}$ = 0 kpc, τ) (Figure 13). If the projected metallicity at the galactic center did not continuously increase, then assuming so would imply non-physical R_birth distributions, with too many stars needing to migrate inwards many kpc in a short time.

Our results show that overall, the assumptions required to estimate [Fe/H](R, τ) and stellar R_birth should hold reasonably well for Milky Way and Andromeda-like disks. The most uncertain assumption is the monotonically enriching $\widehat{[\mathrm{Fe}/\mathrm{H}]}(R_{\mathrm{birth}}=0$ kpc, τ), where we find about half of the galaxies in our sample did not follow this trend. However, our work in Section 6.1.1 suggests that this assumption should hold in the Milky Way, as the R_birth distributions of mono-age populations in the solar neighborhood are meaningful for the Milky Way disk. Looking in detail at a Milky Way-like galaxy (strong bar, similar mass, little to no star formation in the inner disk for the most recent ages), we expect to be able to recover R_birth to within 1 kpc for the Milky Way disk. As age and [Fe/H] measurement uncertainties continue to decrease with future surveys, the ability to estimate R_birth will improve. The increasing precision will also reveal any structure that is currently masked due to observational errors (Renaud et al. 2021a), which should help constrain the recovery of the metallicity profile and stellar birth radii.

While there still exists uncertainty in chemical yields and simulation prescription (Buck et al. 2021), the [Fe/H] gradients in TNG50 are able to reproduce observed MANGA gradients (Lian et al. 2023). Despite having different subgrid physics and implementation, other simulations such as VINTERGATAN (Agertz et al. 2021; Renaud et al. 2021a,b) and NIHAO-UHD (Buck 2020; Buck et al. 2020) are also able to reproduce well global chemical abundance trends such as the chemical bimodality. Therefore, we believe our findings are applicable to observational data and not strongly affected by specific chemical enrichment (not including the dependence on bar strength discovered in this work). Other works have shown that the metallicity gradient can be recovered from the range in [Fe/H] for a few Milky Way-like galaxies using 2 different simulation suites (NIHAO-UHD and HESTIA; Lu et al. 2022a) and additionally on Large Magellanic Cloud-like galaxies using NIHAO-UHD (Lu et al. 2024), thus showing the metallicity gradient-range correlation is not an artifact of TNG50’s prescription. As discussed in Section 5.1, the metallicity profile shows a spike in [Fe/H] near the galactic center. However, as discussed, this drastic increase may be a result of TNG-specific evolution as it does not seem to appear consistently across other simulation suites (e.g. Renaud et al. 2024 and Figure 3 in Lu et al. 2024). The main conclusion of Section 5.1 still holds true though, that regardless of what the metallicity profile looks like at the true galactic center, a projected central metallicity is best to estimate stellar birth radii of disk stars.

Acknowledgments

The authors wish to thank Annalisa Pillepich for her helpful discussions and the Illustris TNG collaboration for making their simulations publicly available, along with an active discussion forum. B.R. and I.M. acknowledge support by the Deutsche Forschungsgemeinschaft under the grant MI 2009/2-1.

References

Appendix A: Galaxies with low correlation between [Fe/H] and R_birth

In Section 4.1 we showed that while most galaxies in the Milky Way and Andromeda analogues catalog of TNG50 (Pillepich et al. 2023) displayed a linear metallicity gradient, there were 17 galaxies that did not have a strong correlation between [Fe/H] and birth radius at a given age. Here, we discuss these galaxies in more detail to illustrate why they do not have a linear gradient, and how they also are non-Milky Way-like. Figure A.1 shows the age–metallicity relation of these 17 galaxies with their subhalo ID number given in the bottom left corner of each panel. The grey colormap corresponds to the density of stellar particles that are currently bound and were born bound to the host galaxy. The colored density corresponds to the density of our stellar disk sample (R < 15 kpc, R_birth < 15 kpc, |z_birth|< 1 kpc, |z|< 1 kpc, eccentricity < 0.5). The dashed lines correspond to merger events with a gas mass ratio > 25% (blue) and 10 − 25% (orange). For 10 of the 17 galaxies, there are minimal stellar disk particles, which span only a few Gyr. 6 of the galaxies (432106, 445626, 452031, 555601, 571633, 588831) exhibit a nonlinear gradient due to a large spike of [Fe/H] in the inner kpc (see Figure 7 for a similar, less severe example). Of the 17 galaxies with poor [Fe/H]–R_birth correlation, Subhlao 571454 appears to be the only one with a non-linear gradient due to a large dispersion in [Fe/H] for a given birth radius at a given time.

Fig. A.1. Age–metallicity relation of the 17 galaxies that were removed from analysis in Section 4 due to the lack of linear metallicity gradient. The grey density map corresponds to the density of stellar particles that are currently bound and were born bound to the host galaxy. The colored density illustrates the density of our stellar disk sample (R < 15 kpc, Rbirth < 15 kpc, |zbirth|< 1 kpc, |z|< 1 kpc, eccentricity < 0.5). The dashed lines show the merger events with a gas mass ratio > 25% (blue) and 10 − 25% (orange). Most of these galaxies are non-Milky Way-like due to too few disk stellar particles. The others are non-Milky Way-like due to a large spike in [Fe/H] in the inner kpc (e.g. subhalo 555601) causing a non-linear metallicity gradient (see also Figure 7).

Appendix B: Extra Figures

Figure B.1shows the relative metallicity gradient for each galaxy as a function of age. We find ∼10 − 15 galaxies show an initial steepening in their gradient before weakening over time, similar as found in the Milky Way (Lu et al. 2022a; Ratcliffe et al. 2023) and in external galaxies.

Figure B.2 shows the relative metallicity gradient versus the range in metallicity as a function of time. The correlation values are used in Section 4.2. Many galaxies have a lower correlation value due to one or two outlier points (e.g. the first panel of row 11).

Figure B.3 illustrates the age–metallicity relation for stellar disk particles (grey density map) and stellar disk particles with R_birth < 1 kpc (blue density) in each galaxy. Many of these galaxies show star formation quenched in their central kpc before redshift 0. Figure B.4 also shows the age–metallicity relation for stellar disk particles in each galaxy, along with the projected central metallicty evolution (blue line). The projected central metallicity is the anchor point for the rest of the disk metallicity profile (Equation 1) and is shown to be more representative of the disk than the inner few kpc (Section 5).

Figure B.5 illustrates how changing the value of the steepest gradient in recovering ∇[Fe/H](τ) from $\mathrm{Range}\stackrel{\sim }{[\mathrm{Fe}/\mathrm{H}]}(\mathrm{age})$ affects the R_birth distributions of mono-age populations in the solar neighborhood. This example is using the incorrect red line from Figure 13 as the central metallicity evolution, and shows how even changing the value of the steepest gradient does not recover the proper distributions as given in the second panel of Figure 13.

Fig. B.1. Evolution of ∇[Fe/H](τ) over time for each galaxy with a linear gradient ([Fe/H]–Rbirth correlation below -0.85) for at least the most recent 4 Gyr of star formation. We find that nearly all galaxies show a weakening in their [Fe/H] gradients over time. About 20 galaxies show a steepening in their metallicity gradient at older ages before flattening, consistent with what is now predicted for the Milky Way.

Fig. B.2. Relationship between scatter in [Fe/H] and [Fe/H] radial birth gradient. Each point represents a different age bin, colored from youngest (blue) to oldest (yellow). The grey line is the line of best-fit. The correlation between ∇[Fe/H](τ) and $\mathrm{Range}\stackrel{\sim }{[\mathrm{Fe}/\mathrm{H}]}(\mathrm{age})$ is given in the bottom left of each panel. The red font indicates the galaxies analysed in Figure 5.

Fig. B.3. Age–metallicity relation of disk stellar particles for each galaxy in our sample. The grey density refers to the disk sample, while the blue density corresponds to the density of the disk sample with Rbirth < 1 kpc. Most galaxies show a halt in their star formation in their central regions.

Fig. B.4. Density of the age–metallicity relation for the disk stellar sample of each galaxy. The dark blue line corresponds to the projected central metallicity that is used as the anchor point to determine the metallicity profile for the genuine disk (see Section 5). About half the sample shows a dilution in the projected central metallicity, while the other half shows a predominantly [Fe/H]-enriching center. If no blue line is shown, then there were not enough stars to determine the projection.

Fig. B.5. Assuming different gradient strengths (in dex/kpc) changes the shape of the mono-age Rbirth populations for disk stellar particles currently near 8 kpc in subhalo 560751. However, varying the maximum strength of the gradient is never able to recover the correct trends as found in the second from the left panel of Figure 13.

All Figures

Fig. 1. Median absolute correlation between [Fe/H] and Rbirth for galaxies in the TNG50 Milky Way and Andromeda catalog across lookback time grouped by bar strength, with the shaded areas showing ±1 standard deviation about the median. The correlation between [Fe/H] and Rbirth was calculated for age bins of 1 Gyr. To accommodate for outliers and higher feedback, the correlation was calculated for the running median, using Rbirth bins of 0.25 kpc. We do not find that bar strength has an impact on the linearity of the average birth metallicity gradient. The large variability in the correlation between [Fe/H] and Rbirth for stronger barred galaxies at smaller lookback times is predominantly due to weak metallicity gradients, even though the gradient itself is linear.

In the text

Fig. 2. Example (subhalo 577372) illustrating the recovery of ∇[Fe/H](τ) from $\mathrm{Range}\stackrel{\sim }{[\mathrm{Fe}/\mathrm{H}]}(\mathrm{age})$. Left: Density distribution of the [Fe/H]–Rbirth plane for stars born 2 − 3 Gyr ago. The correlation between the [Fe/H] running median and Rbirth is given in the top right corner. The fitted true birth gradient is shown as the black line, and the maximum and minimum [Fe/H] used to determine the range are shown as the dashed red lines. Right: ∇[Fe/H](τ) versus $\mathrm{Range}\stackrel{\sim }{[\mathrm{Fe}/\mathrm{H}]}(\mathrm{age})$ for the same galaxy, with each point representing a different age bin (width 1 Gyr). $\mathrm{Range}\stackrel{\sim }{[\mathrm{Fe}/\mathrm{H}]}(\mathrm{age})$ is determined by looking at the range in [Fe/H] across all age bins, and then normalized so the minimum range is given a value of 0, and the maximum range is given a value of 1. The correlation between $\mathrm{Range}\stackrel{\sim }{[\mathrm{Fe}/\mathrm{H}]}(\mathrm{age})$ and ∇[Fe/H](τ) for this galaxy is −0.87, indicating the evolution of the gradient can be recovered well from $\mathrm{Range}\stackrel{\sim }{[\mathrm{Fe}/\mathrm{H}]}(\mathrm{age})$.

In the text

Fig. 3. Density distribution of the correlation between $\mathrm{Range}\stackrel{\sim }{[\mathrm{Fe}/\mathrm{H}]}(\mathrm{age})$ and ∇[Fe/H](τ) for 179 galaxies (black). These galaxies are chosen such that they have a linear gradient for (at least) their most recent 4 Gyr of star formation. Most galaxies exhibit a reasonable ability to have their ∇[Fe/H](τ) evolution recovered from $\mathrm{Range}\stackrel{\sim }{[\mathrm{Fe}/\mathrm{H}]}(\mathrm{age})$ (correlation ∼ − 1). This ability has a correlation with bar strength (colored lines), where stronger barred galaxies have a better recovery rate than weaker ones.

In the text

Fig. 4. Left: Star-forming region for 2 example galaxies (subhalos 506720 and 539333). The shaded region at a specific age indicates the radial range of stars forming during that time. The star-forming region of some galaxies (such as 539333) can vary in width of many kpc, while other galaxies (such as 506720) form stars in a much more consistent range. Right: Correlation between the scatter in [Fe/H] and the metallicity birth gradient as a function of how much the radial range of the star formation (RRSF) varies over time (defined as the difference between the maximum and minimum range stars are formed in the disk over time), colored by the galaxy’s stellar mass within 30 kpc. A value of 0 indicates that the RRSF was constant over time (albeit potentially covering different radii such as subhalo 506720 in the left panel), while a larger value indicates that the galaxy formed stars over a few kpc some Gyr, and over a much larger region other Gyr (such as subhalo 539333 in the left panel). Overall, the more the range of the star-forming region varies over time, the less likely ∇[Fe/H](τ) can be recovered from $\mathrm{Range}\stackrel{\sim }{[\mathrm{Fe}/\mathrm{H}]}(\mathrm{age})$. There is also a gradient with stellar mass, with higher mass galaxies showing a better ability to have their metallicity gradients recovered from $\mathrm{Range}\stackrel{\sim }{[\mathrm{Fe}/\mathrm{H}]}(\mathrm{age})$. The 4 highlighted galaxies in red boxes are examined in detail in Figure 5.

In the text

Fig. 5. Radial range of star-forming region (RRSF) versus the relative strength of the metallicity gradient for different ages, for the four highlighted galaxies in Figure 4. The correlation between ∇[Fe/H](τ) and $\mathrm{Range}\stackrel{\sim }{[\mathrm{Fe}/\mathrm{H}]}(\mathrm{age})$ is given in the bottom left of each panel, below the subhalo ID. The top panels depict situations in which the evolution of the metallicity gradient can be recovered from the scatter in [Fe/H], while the bottom row illustrates two examples of galaxies where the gradient cannot be recovered in this fashion due to the gradient being steepest during a time of a lower star-forming region.

In the text

Fig. 6. Top: Growth of the regions where stars form in the galactic disk over time, separated by bar strength at redshift 0. Weaker-barred galaxies start with stars forming on less than 3 kpc in the center of the disk, and then grow until stars are forming almost throughout the entirety of the disk. Stronger-barred galaxies show only minimal growth in their radial range of star-forming region (RRSF) Middle: Upper radial limit of where stars are forming throughout the galaxy at a given moment. Each group of galaxies shows the expected inside-out growth. Bottom: Lower radial limit of where stars are forming in the galaxy over time. Weaker-barred galaxies tend to form stars consistently in their inner regions, whereas stronger-barred galaxies stop forming stars in the inner few kpc. In each panel, the lines correspond to the median of each group and the shaded region represents the 25−75% quantiles. To best illustrate the growth of the disk over time, we use age bins of 0.2 Gyr.

In the text

Fig. 7. Illustrating the recovery of [Fe/H](Rbirth = 0 kpc) for subhalo 553837 for stellar particles 7−8 Gyr ago. The projected central metallicity ($\widehat{[\mathrm{Fe}/\mathrm{H}]}(R_{\mathrm{birth}}$ = 0 kpc); black star) was estimated from the metallicity at a given radius (8.5 kpc; black square) and the metallicity gradient at the current time. The red circle represents the true central metallicity [Fe/H](Rbirth = 0 kpc), taken as the median [Fe/H](Rbirth < 1 kpc). For this scenario, the central metallicity cannot be estimated accurately assuming a linear gradient due to the spike in [Fe/H] in the inner kpc. This galaxy has a weak bar (bar strength of 0.15) and a mass of 3.5 × 1010 M⊙.

In the text

Fig. 8. Density distribution of the difference between inferred and true [Fe/H] at galactic center separated by top: lookback time τ and bottom: bar strength. The projected [Fe/H](R = 0 kpc, τ) is estimated from the [Fe/H] value at a given radius and the current [Fe/H] gradient for a given age of width 1 Gyr (see Figure 7). As a whole, the value of [Fe/H](R = 0 kpc, τ) is underestimated by ∼0.2 dex, with larger differences for weaker barred galaxies and larger lookback times. The bias towards negative values (i.e. [Fe/H](R = 0 kpc, τ) is underestimated) is predominantly due to the non-linear spike in [Fe/H] in the inner kpc (see Figure 7 for an example).

In the text

Fig. 9. Schematic illustrating the distribution of galaxies and how their left: [Fe/H](R = 0 kpc, τ) and right: $\widehat{[\mathrm{Fe}/\mathrm{H}]}(R_{\mathrm{birth}}$ = 0 kpc) evolves over time. The evolution of the true central metallicity (left) was estimated by tracing how [Fe/H] in the central kpc of the galactic disk sample enriches over time. Of the 183 galaxies that have enough disk stars born within the inner kpc, only 15 of them are found to have their centers predominantly continuously enrich in [Fe/H] over cosmic time, while 45 of the galaxies are found to have at least one major depletion event in their galactic centers. The remaining 123 galaxies had star formation quench before redshift 0 in their center. An example of each of these scenarios is given, where the density of the age–metallicity relation for our defined disk sample (grey) and stellar particles with Rbirth < 1 kpc (blue) is shown. The evolution of the projected central metallicity (right) was determined by tracing the projected [Fe/H] at Rbirth = 0 kpc using the metallicity at a given radius and the current metallicity gradient (see Figure 7 for an example). $\widehat{[\mathrm{Fe}/\mathrm{H}]}(R_{\mathrm{birth}}$ = 0 kpc) traces the central metallicity assuming a linear metallicity gradient, and acts as an anchor point for the disk’s gradient, regardless if the true central metallicity spikes in the galactic center. Of the 178 galaxies with enough disk stars to appropriately trace the time evolution of $\widehat{[\mathrm{Fe}/\mathrm{H}]}(R_{\mathrm{birth}}$ = 0 kpc), we find 81 galaxies show a continuous enrichment at Rbirth = 0 kpc, while 97 have at least one major depletion event. Examples of these two scenarios are given next to the respective outcome, with the blue line indicating the projected central metallicity, estimated using age bins of 0.5 Gyr. The central metallicity of the Milky Way is assumed to continuously increase in current methods, with potentially only minor deviations from this assumption (see Section 6.1.1).

In the text

Fig. 10. Distribution of times between the depletion event in $\widehat{[\mathrm{Fe}/\mathrm{H}]}(R_{\mathrm{birth}}$ = 0 kpc, τ) and the last major (blue; gas mass ratio > 25%) and minor/major (orange; gas mass ratio > 10%) mergers. A depletion event is defined as the projected central metallicity decreasing by at least 0.05 dex in 0.5 Gyr, smoothed over 2.5 Gyr. The majority of the time that the projected metallicity in the galactic central region decreases, the galaxy experienced a merger of at least 10% gas mass fraction within the previous 2 Gyr. This suggests that larger mergers potentially affect the evolution of the metallicity of the galactic center.

In the text

Fig. 11. Density of the change in $\widehat{[\mathrm{Fe}/\mathrm{H}]}(R_{\mathrm{birth}}$ = 0 kpc, τ) as a function of time after a (left) massive or (right) minor merger for all galaxies after all massive/minor mergers, until either another massive/minor merger occurs or no more stars form in the central 1 kpc. A positive change in [Fe/H] indicates that the projected metallicity in the central 1 kpc increased, while a negative value indicates that the projected [Fe/H] in the central kpc decreased. The projected metallicity at the galactic center does not necessarily decrease immediately after a merger, and a dilution in $\widehat{[\mathrm{Fe}/\mathrm{H}]}(R_{\mathrm{birth}}$ = 0 kpc, τ) can take up to a few Gyr before diluting the projected galactic center in [Fe/H]. Even if $\widehat{[\mathrm{Fe}/\mathrm{H}]}(R_{\mathrm{birth}}$ = 0 kpc, τ) does not dilute, there is a general trend that [Fe/H] enrichment slows down.

In the text

Fig. 12. Two examples of the evolution of $\widehat{[\mathrm{Fe}/\mathrm{H}]}(R_{\mathrm{birth}}$ = 0 kpc, τ) after a merger with pre-enriched gas. Due to the metal-rich spike in the galactic center (Figure 7), we remove the inner kpc for visualization purposes. The subhalo in the top row (subhalo 554523) had a major merger at lookback time ∼9 Gyr ago which brought in more metal-rich gas and increased the metallicity of the entire disk. This allowed for the projected central metallicity (marked with the red “x”) to also increase. The subhalo in the bottom row (subhalo 514829) experienced a metal-rich merger ∼8 Gyr ago. However, instead of increasing the metallicity throughout the entire disk like in the top row, here only the outer disk had an increase in [Fe/H]. This thus weakened the metallicity gradient, and in turn lowered the projected central metallicity. These examples show that even something as trivial as the metallicity of the merger compared to the main galaxy cannot begin to disentangle when $\widehat{[\mathrm{Fe}/\mathrm{H}]}(R_{\mathrm{birth}}$ = 0 kpc, τ) increases or decreases after a merger.

In the text

Fig. 13. Ability to recover Rbirth for subhalo 560751 with age and [Fe/H] uncertainties added (15% and 0.05 dex respectively). Left: The age–metallicity relation for this subhalo with the true evolution of $\widehat{[\mathrm{Fe}/\mathrm{H}]}(R_{\mathrm{birth}}$ = 0 kpc, τ) (black line) and estimated evolution assuming the projected galactic center is always increasing in [Fe/H] (red line). The second, third, and right panels show the Rbirth distributions for mono-age populations of solar neighborhood (7.75 ≤ R ≤ 8.25) stellar particles for the true Rbirth, Rbirth inferred using the true $\widehat{[\mathrm{Fe}/\mathrm{H}]}(R_{\mathrm{birth}}$ = 0 kpc, τ) (InferreddecRbirth), and Rbirth inferred using the monotonically increasing $\widehat{[\mathrm{Fe}/\mathrm{H}]}(R_{\mathrm{birth}}$ = 0 kpc, τ) (InferredincRbirth) respectively. Inappropriately assuming a monotonically increasing metallicity abundance in the galactic center provides Rbirth distributions that are non-physical and do not match the true distributions, whereas using the true $\widehat{[\mathrm{Fe}/\mathrm{H}]}(R_{\mathrm{birth}}$ = 0 kpc, τ) evolution provides a more accurate trend in how Rbirth decreases with increasing age. Examining the Rbirth distributions of mono-age populations in the solar neighborhood is thus a necessity in constraining $\widehat{[\mathrm{Fe}/\mathrm{H}]}(R_{\mathrm{birth}}$ = 0 kpc, τ) evolution.

In the text

Fig. 14. Inferred Rbirth versus true Rbirth for subhalo 560751 using the method from Lu et al. (2022a) with age and [Fe/H] uncertainties added of 15% and 0.05 dex respectively. Rbirth are estimated using the true evolution of [Fe/H](R = 0 kpc, τ) (black line in Fig. 13) and ∇[Fe/H](τ) recovered from $\mathrm{Range}\stackrel{\sim }{[\mathrm{Fe}/\mathrm{H}]}(\mathrm{age})$. The difference between the left and right panels is that stellar [Fe/H] values in the TNG50 galaxy were adapted to better match the dispersion about the metallicity gradient as seen in the Milky Way’s H II regions in the right panel. This illustrates the high probability of success in estimating Rbirth in the Milky Way disk.

In the text

Fig. 15. Left column: Evolution of mono-Rbirth populations in the age–metallicity relation for a TNG50 Milky Way and Andromeda-like galaxy (subhalo 402555) whose mono-Rbirth populations exhibit “wiggling” behavior, where each mono-Rbirth population did not experience a monotonic enrichment in [Fe/H]. Right column: [Fe/H] abundance of mono-Rbirth populations relative to [Fe/H](R < 2 kpc) across age. Looking at the [Fe/H] abundance relative to the inner few kpc at a given age allows the structure of the wiggles to become more clear. Most of the galaxies that exhibit “wiggling” structure in this plane follow similar trends to this galaxy – each Rbirth sees a depletion in [Fe/H] at the same time as also seen using Rbirth methods (see Fig. 2 in Ratcliffe et al. 2023).

In the text

Fig. A.1. Age–metallicity relation of the 17 galaxies that were removed from analysis in Section 4 due to the lack of linear metallicity gradient. The grey density map corresponds to the density of stellar particles that are currently bound and were born bound to the host galaxy. The colored density illustrates the density of our stellar disk sample (R < 15 kpc, Rbirth < 15 kpc, |zbirth|< 1 kpc, |z|< 1 kpc, eccentricity < 0.5). The dashed lines show the merger events with a gas mass ratio > 25% (blue) and 10 − 25% (orange). Most of these galaxies are non-Milky Way-like due to too few disk stellar particles. The others are non-Milky Way-like due to a large spike in [Fe/H] in the inner kpc (e.g. subhalo 555601) causing a non-linear metallicity gradient (see also Figure 7).

In the text

Fig. B.1. Evolution of ∇[Fe/H](τ) over time for each galaxy with a linear gradient ([Fe/H]–Rbirth correlation below -0.85) for at least the most recent 4 Gyr of star formation. We find that nearly all galaxies show a weakening in their [Fe/H] gradients over time. About 20 galaxies show a steepening in their metallicity gradient at older ages before flattening, consistent with what is now predicted for the Milky Way.

In the text

Fig. B.2. Relationship between scatter in [Fe/H] and [Fe/H] radial birth gradient. Each point represents a different age bin, colored from youngest (blue) to oldest (yellow). The grey line is the line of best-fit. The correlation between ∇[Fe/H](τ) and $\mathrm{Range}\stackrel{\sim }{[\mathrm{Fe}/\mathrm{H}]}(\mathrm{age})$ is given in the bottom left of each panel. The red font indicates the galaxies analysed in Figure 5.

In the text

	Fig. B.3. Age–metallicity relation of disk stellar particles for each galaxy in our sample. The grey density refers to the disk sample, while the blue density corresponds to the density of the disk sample with Rbirth < 1 kpc. Most galaxies show a halt in their star formation in their central regions.
In the text

Fig. B.4. Density of the age–metallicity relation for the disk stellar sample of each galaxy. The dark blue line corresponds to the projected central metallicity that is used as the anchor point to determine the metallicity profile for the genuine disk (see Section 5). About half the sample shows a dilution in the projected central metallicity, while the other half shows a predominantly [Fe/H]-enriching center. If no blue line is shown, then there were not enough stars to determine the projection.

In the text

	Fig. B.5. Assuming different gradient strengths (in dex/kpc) changes the shape of the mono-age Rbirth populations for disk stellar particles currently near 8 kpc in subhalo 560751. However, varying the maximum strength of the gradient is never able to recover the correct trends as found in the second from the left panel of Figure 13.
In the text