A Tale of Two Disks: Mapping the Milky Way with the Final Data Release of APOGEE

The Apache Point Observatory Galactic Evolution Experiment (APOGEE; Majewski et al., 2017) is a high-resolution ($R\sim 22,500$) near-infrared ($1.51-1.70$ $\mu$m) spectroscopic survey containing observations of 657,135 unique stars released as part of the SDSS-IV survey (Blanton et al., 2017). The spectra were obtained using the APOGEE spectrograph (Wilson et al., 2019) mounted on the $2.5$m SDSS telescope (Gunn et al., 2006) at Apache Point Observatory to observe the Northern Hemisphere (APOGEE-N), and expanded to include a second APOGEE spectrograph on the $2.5$ m Irénée du Pont telescope (Bowen & Vaughan, 1973) at Las Campanas Observatory to observe the Southern Hemisphere (APOGEE-S). The final version of the APOGEE catalog was published in December 2021 as part of the 17th data release of the Sloan Digital Sky Survey (DR17; Abdurro’uf et al., 2022) and is available publicly online through the SDSS Science Archive Server and Catalog Archive Server²²2Data Access Instructions: https://www.sdss.org/dr17/irspec/spectro_data/.

The APOGEE data reduction pipeline is described in Nidever et al. (2015). Stellar parameters and chemical abundances in APOGEE were derived within the APOGEE Stellar Parameters and Chemical Abundances Pipeline (ASPCAP; Holtzman et al., 2015; García Pérez et al., 2016; Holtzman et al., 2018; Jönsson et al., 2020, J.A. Holtzman et al. 2022 in prep.). ASPCAP derives stellar atmospheric parameters, radial velocities, and as many as 20 individual elemental abundances for each APOGEE spectrum by comparing each to a multi-dimensional grid of theoretical model spectra (Mészáros et al., 2012; Zamora et al., 2015) and corresponding line lists (Shetrone et al., 2015; Smith et al., 2021), employing a $\chi^{2}$ minimization routine with the code FERRE (Allende Prieto et al., 2006) to derive the best-fit parameters for each spectrum. We highlight that several elements (notably [Mg/Fe]) were updated in DR17 to include non-LTE effects in the stellar atmosphere. ASPCAP reports typical accuracy in metallicity measurements within 0.01 dex (Jönsson et al., 2018). In this study, we adopt the calibrated values for surface gravity ($\log g$), metallicity ([Fe/H]), and $\alpha$-element abundances ([Mg/Fe]) from ASPCAP. We adopt [Mg/Fe] for our $\alpha$-element abundance instead of the “total” [$\alpha$/M], because [Mg/Fe] is the most precisely measured abundance by ASPCAP, and this element ratio has been traditionally used to define the boundary between the chemical thin and thick disk.

Several cuts were made to the full APOGEE catalog to refine our sample. First, only stars defined as APOGEE main survey targets (also sometimes called the “main red giant sample”) were selected using the EXTRATARG flag. This removes any duplicate entries, as well as any ancillary science or other survey stars that were targeted for observation for a specific purpose (e.g., satellite or dwarf galaxy targets, star cluster member candidates, Kepler Objects of Interest). The main survey targets were randomly selected for observation from the 2MASS catalog, based on their $(J-K)$ color and $H$-band apparent magnitude. For more information on the targeting strategies of APOGEE, see Zasowski et al. (2013, 2017), Beaton et al. (2021), and Santana et al. (2021).

Stars with noisy spectra ($S/N<50$) or unreliable parameter estimates from ASPCAP were removed from our sample using the SN_BAD and STAR_BAD ASPCAP bits respectively. The latter is triggered when the derived parameters for a star are designated a bad fit by its high $\chi^{2}$ value, when the derived temperature does not match the star’s observed color, when any individual stellar parameter measurement is flagged as bad, or when the derived parameters lie on an edge of the synthetic spectral grid.

The sample is further restricted to stars with surface gravity values between $1\leq\log g\leq 2$. Limiting to a small range in $\log g$ minimizes potential systematic uncertainties in abundance measurements, which tend to present as a function across T_eff and $\log g$ in APOGEE (e.g., Jönsson et al., 2018; Eilers et al., 2021). The higher luminosity of these giants helps probe larger distances, allowing for a wide range of positions to be sampled across the Galactic disk in our study. Fainter stars may be better sampled closer to the Sun, but to keep our sample consistent across all distances, we apply this $\log g$ cut to ensure the trends we are documenting are not attributed to any systematic bias. Eilers et al. (2021) presents an empirical correction for these systematics for those interested; this would mainly be a concern if the expected distribution of $\log g$ in observations varies significantly with distance, which we do not expect in our sample. The lower $\log g$ limit is also imposed by the availability of asteroseismic data, as no age estimates are available for stars with $\log g<1$ in distmass.

Figure 1 shows the T_eff-$\log g$ distribution of the sample after these refinements. The final number of stars in our RGB sample is 66,496. Due to the particulars of APOGEE field selection, there are more stars above the disk ($Z\geq 0;N=$ 38,031) than below ($Z\leq 0;N=$ 28,465), and more observations towards the Galactic center ($R\leq R_{\odot};N=$ 36,317) than outward ($R\geq R_{\odot};N=$ 30,179). The spatial distribution of our sample is shown in Figure 2. The stellar distance estimates used for this Figure (and the remaining of the paper) are described in Section 2.4.

In this work, we make no correction for the selection biases within the APOGEE survey. Stars close to the solar neighborhood will be over-represented in our sample. As shown in Figure 2, as distance from the Sun increases, the number of stars available in the APOGEE sample decreases.³³3This explanation is an oversimplification of the APOGEE selection function. The actual selection function depends on more than just distance from the Sun, as targeting strategies may vary between fields, observing time availability and instrument specifications differed between the North and South, and the non-homogeneous dust distribution in the Milky Way plays a major role in what can be observed. There are certain limitations that this selection function imposes on this work and similar studies. Specifically, results should not be averaged over a large spatial range, as the relative number of observed stars will clearly weight the average towards the solar neighborhood. Additionally, nothing can be inferred from the relative number of stars between locations, or about the intrinsic density profile of stars in the disk, because the former is heavily influenced by the selection function. That said, the effect of the selection function should be negligible when confined within a small spatial zone and $\log g$ limit in the Galaxy, such that general abundance trends and normalized number distributions should be consistent even without correcting for the APOGEE selection function (e.g., Hayden et al., 2015, Appendix A). In this work, we consistently bin stars in different ranges of $R$ and $Z$ to avoid this bias and explore how chemical and age trends vary across the disk. A more complete prescription on how to account for the effects of the selection function in APOGEE has been published in Bovy et al. (2012b, 2016), Mackereth et al. (2017) for previous data releases, and Imig et al. (in prep.) for DR17. Imig et al. (in prep.) will present the density distribution of mono-age mono-abundance stellar populations in APOGEE DR17 after correcting for the selection function.

Figure 3 shows the [Mg/Fe]-[Fe/H] plane for our full sample, where the bimodal distribution in [Mg/Fe] is obvious (e.g., Fuhrmann, 1998; Bensby et al., 2005; Reddy et al., 2006; Lee et al., 2011; Kordopatis et al., 2015a; Hayden et al., 2015; Katz et al., 2021). Although there is some debate on whether the two sequences are truly distinct (see Introduction), we use this figure to define two further subsamples in our data to investigate this question later. Vincenzo et al. (2021a) demonstrate that the distribution in [Mg/Fe] at fixed [Fe/H] is genuinely bimodal when considering the full disk population at near-solar radii, after accounting for the APOGEE selection function. We define the $\alpha$-poor “thin disk” sequence and the $\alpha$-rich “thick disk” sequence by splitting the full sample into two groups defined by a line in [Mg/Fe]-[Fe/H] space, shown in Figure 3. We adopt a similar limit as Weinberg et al. (2019, 2022) parameterized by the equation:

Our equation differs from Weinberg et al. (2019, 2022) by a small offset of [Mg/Fe] $=+0.03$ dex, correcting for abundance calibrations.

This separation in $\alpha$-element abundances is shown in Figure 3. A conservative buffer zone within $\pm 0.025$ dex of the line is excluded to remove potential overlap between the two sequences. This value is larger than the typical uncertainties of APOGEE abundance measurements, shown as the $\pm 1\sigma$ in the error bars across the bottom of the plot.

Accurately mapping the Milky Way in three dimensions requires knowing precise distances to every star in our sample. Galactocentric positions were calculated for each star using the right ascension (RA) and declination (DEC) from APOGEE observations and distance estimates from the APOGEE distmass value added catalog (A. Stone-Martinez et al. 2023, submitted). The distmass distances were obtained through a neural network that was trained to estimate a star’s luminosity based on its ASPCAP parameters, using Gaia and cluster distances to provide the training labels. Distance estimates from the distmass catalog are typically precise within 10%. For the purpose of calculating Galactocentric coordinates, we define the reference location of the Sun to be $R_{\odot}=8.3$ kpc with a height of $z_{\odot}=0.027$ kpc above the plane (Bland-Hawthorn & Gerhard, 2016).

For evolved red giant stars, carbon and nitrogen abundances provide mass information because of the mass-dependence of stellar mixing (e.g., Iben, 1965; Salaris & Cassisi, 2005), allowing the determination of stellar masses for stars without asteroseismology (e.g., Martig et al., 2016; Vincenzo et al., 2021b). This fundamental property is used to derive stellar age estimates in the distmass catalog, wherein ages are derived by training a second neural network on the ASPCAP parameters of stars with asteroseismology masses from the APOGEE-Kepler overlap survey (APOKASC; Pinsonneault et al., 2018, APOKASC3: Pinsonneault et al. 2023 in prep.). The neural network learns the relations between the ASPCAP parameters and asteroseismic masses for stars from APOKASC, then it predicts the masses for all giant stars from DR17. Knowing the masses for evolved stars, ages can be derived through stellar evolution theory which predicts a star’s location on an isochrone. For the distmass catalog, isochrones from Choi et al. (2016) were adopted to make this conversion from derived mass to stellar age. The isochrones cover a range of ages ($5.0\leq\mathrm{log}(\mathrm{age})\leq 10.3$), metallicities ($-2.0\leq[{\rm{Z}}/{\rm{H}}]\leq 0.5$), and masses ($0.1\leq M/{M}_{\odot}\leq 300$).

The reported uncertainties on the age estimates in distmass are shown as a histogram in Figure 4 for our sample; the median lower uncertainty is $0.14$ $\log_{10}$(age) and the median upper uncertainty is $0.16$ $\log_{10}$(age). The uncertainties in stellar age are propagated from the uncertainties in stellar mass, which are predicted from the spread in mass from the neural network training set. The mass and age uncertainties have no strong dependence on metallicity or $\log g$ (A. Stone-Martinez et al. 2023, submitted).

To determine if these uncertainties are realistic, A. Stone-Martinez et al. (2023, submitted) performs an additional evaluation of the age estimates by comparing to previous literature. Compared to a small sample of cluster members with independent age estimates from main sequence turnoff fitting, they find that the distmass ages are accurate within $\pm 1\sigma=0.16$ in $\log_{10}$(age) across 12 different star clusters with ages $9.2\leq\log_{10}$(age) $\leq 9.7$ Gyr. Compared to a larger sample of field stars from the astroNN catalog (Leung & Bovy, 2018; Mackereth et al., 2019), the distmass ages show a typical spread of $\pm 1\sigma=0.11$ in $\log_{10}$(age), although these ages are derived with a similar methodology. Both of these evaluations are consistent within the reported age uncertainties.

Precise stellar ages remain difficult to measure robustly for large samples of stars. Neural network derived ages like distmass are heavily dependent on their training set values, and any uncertainties in the training set labels will influence the neural network model. The distmass catalog trains on stellar masses derived from asteroseismology, which have their own systematics and different groups have derived different results. Appendix A tests some of our results using different training sets and age catalogs, and motivates our choice of adopted ages. Even within distmass, there are six provided stellar age estimates for every star from neural networks trained on different asteroseismic mass estimates from three different research groups in APOKASC 3 (Pinsonneault et al. 2023 in prep.), and a set of “corrected” and “uncorrected” masses for each. The corrections are motivated by Gaia data, using parallaxes to calibrate the stellar radius derived from asteroseismology (e.g., Zinn et al., 2019). These corrections may be less reliable at lower $\log g$ and produce results that are less consistent across different bins in $\log g$. For the remainder of this paper, we use the distmass results trained on the uncorrected SS ages from APOKASC 3 (Pinsonneault et al. 2023 in prep.), corresponding to the column named “AGE_UNCOR_SS” in the distmass catalog.

Selecting different age catalogs does make some difference in our results, particularly among the oldest stars as shown in Appendix A. Because this quantitative aspect can change considerably, we caution the reader against drawing strong conclusions from any of our age-related results without thoroughly understanding the related caveats outlined in Appendix A.

An additional quality flag from the distmass catalog is used to refine the sample when using the stellar age estimates. Namely, we remove stars that have bit 2 set, indicating stellar parameters T_eff, [Fe/H], and $\log g$ lie outside of the range covered by the APOKASC training set; this removes stars with potentially unreliable mass (and therefore age) estimates. For anything involving ages, the full RGB sample is additionally restricted to 57,756 stars with this cut. Notably, all stars with [Fe/H] $\leq-0.7$ are excluded by this criterion. Metal-poor stars have extra mixing that was not learned by the neural network because there were no metal poor stars in the training set. Because metal-poor stars tend to be older, this means that for age related figures, the oldest stars ($\rm{age}>10^{10}$ years) may not be well-represented in our sample, particularly at large radii in the Galaxy; the potential effects of this, and other age-related caveats, are explored more in Appendix A.

Maps of the Galactic disk as a function of [Fe/H], [Mg/Fe], and stellar age are shown in Figure 5 in face-on ($X-Y$ plane; top row) and edge-on ($R-Z$ plane; bottom row) perspectives. The median value of each parameter is calculated for different spatial bins sized $\Delta X=\Delta Y=0.5$ kpc, and shown as the respective color on the figure. For the edge-on perspective, the sign of the $X$ coordinate is applied to the $R$ coordinate, to better highlight the spatial coverage of the observations on the opposite side of the Galaxy.

In median metallicity (left column), clear radial and vertical metallicity gradients are visible in the disk, with higher average metallicities near the Galactic center that decline towards outer radii. In median $\alpha$-element abundances (middle column), the bimodality in the disk shows low-$\alpha$ stars congregating in the “thin disk” near the Galactic midplane, and high-$\alpha$ stars populating the “thick disk” at higher $Z$ locations. At larger radii, the low-$\alpha$ stars extend farther above and below the plane. The high-$\alpha$ stars are more centrally concentrated. The right column, colored by median stellar age, contains fewer stars due to the additional cuts described in Section 2.4 when dealing with ages from the distmass catalog. Once again, radial and vertical gradients appear in these maps, as well as younger stars extending farther above the plane in the outer Galaxy. The innermost structural features of the Galaxy, such as the bar (noted by the ellipse) and bulge stand out as metal-rich, $\alpha$-poor, and older-aged than stars at similar radii but different azimuthal angles, consistent with previous studies of the central regions of the Galaxy (Wegg et al., 2019; Zasowski et al., 2019; Hasselquist et al., 2020; Eilers et al., 2021; Lian et al., 2021; Queiroz et al., 2021) for this metallicity range.

Dividing the maps into vertical bins reveals more nuanced structure; Figure 6 depicts face-on metallicity maps divided by height above the Galactic plane, from closest to the Galactic plane (top panel; $|Z|\leq 0.5$ kpc) to farthest away (bottom panel; $1\leq|Z|\leq 2$ kpc). The metallicity gradient is strongest close to the Galactic plane, with locations near the Galactic center showing a higher median metallicity than those at larger radii, as expected. Farther above the mid plane, the trend becomes less apparent, with almost no obvious gradient present when $|Z|\geq 1$ kpc, and the stellar populations showing a lower median metallicity overall. The middle panel ($0.5\leq|Z|\leq 1$ kpc) shows a peculiar trend where the median metallicity actually increases from $R=0$ until $R\sim 7$ kpc, and then decreases with a shallow metallicity gradient. This build-up of metal-rich stars in the center of the Galaxy is possibly the signature of the bulge.

The age distribution of the Galactic disk is shown in Figure 7. Again, the age gradient is strong close to the Galactic plane, with older stars more common near the center and younger stars dominating in the outer Galaxy. Unlike in metallicity, there is still a radial age gradient above the Galactic plane ($|Z|>1$ kpc), although in general the stars found above the plane are older than the stars found in the plane.

Figure 8 shows the Galactic distribution of stars in the [Mg/Fe]-[Fe/H] chemistry plane as a function of Galactic position. The different rows are the same vertical bins adopted in previous sections, with the bottom row closest to the Galactic plane ($0<|Z|<0.5$ kpc) and the top row farthest from the plane ($1.0<|Z|<2.0$ kpc). The columns are different radial bins, from closest to the disk center (left column) to farthest out (right column). Each panel shows the [Mg/Fe]-[Fe/H] distribution of stars in its respective spatial zone, colored by stellar point density (in Figure 8) and stellar age (in Figure 9). Our adopted definition of the split between the high- and low-$\alpha$ sequences (equation 1) is plotted in black. For reference, the gray background highlights the distribution of the full sample, indicating the contour within which 90% of the sample is found.

Generally, the low-$\alpha$ sequence is concentrated close to the Galactic plane (bottom row), and the high-$\alpha$ sequence is more prominent outside the plane (top row) for $R<12$ kpc. The location of the high-$\alpha$ sequence does not change based on location in the Galaxy. The low-$\alpha$ sequence is more metal rich near the center of the Galaxy (left column), and moves to more metal-poor with increasing radius (right column). Additionally, at large radii, the low-$\alpha$ sequence extends farther above the plane than it does close to the Galactic center. All of this has been well-documented in previous studies (e.g., Bensby et al., 2005, 2011; Nidever et al., 2014; Kordopatis et al., 2015a; Hayden et al., 2015; Katz et al., 2021; Vincenzo et al., 2021a; Gaia Collaboration et al., 2022).

Figure 9 shows the spatial distribution of [Mg/Fe] coded by stellar age. The high-$\alpha$ sequence is composed of older stars at all radial bins. The low-$\alpha$ sequence includes older stars close to the Galactic center ($R<6$ kpc) and younger stars farther out in radius. At any radius, the lower [Mg/Fe] stars within the low-$\alpha$ sequence have younger ages. Within the low-$\alpha$ sequence, stellar age correlates more with [Mg/Fe] than it does with [Fe/H], indicating that the low-$\alpha$ sequence is likely not a true single sequence.

To aid in the direct comparison between the radial and height bins, Figure 10 shows the contours (top panel) and median (bottom panel) in the [Mg/Fe]-[Fe/H] plane for both the low-$\alpha$ and high-$\alpha$ samples as a function of Galactic radius. In this Figure, the data are restricted to the plane ($|Z|<0.5$ kpc), equivalent to the bottom row in Figures 8 and 9. The contour containing 90% of points in both the low-$\alpha$ and high-$\alpha$ samples is shown in the top panel. The low-$\alpha$ sequence is more metal-rich near the center of the Galaxy, and shifts continuously to lower metallicities and higher-$\alpha$ moving outwards in radius. The high-$\alpha$ sequence’s contour is generally the same shape and position at all radii, although close to the center of the Galaxy, the shape extends further to the metal-poor end. In the bottom panel, the median [Mg/Fe] as a function of metallicity is shown for both samples. As before, the low-$\alpha$ sample shifts more metal-poor with increasing radius. The high-$\alpha$ sequence moves slightly downwards (towards more $\alpha$-poor) with increasing radius. This is also seen in Katz et al. (2021) with APOGEE data, using the mode of the data, although they find a larger shift of $\sim 0.05$ dex between the inner and outer Galaxy, while ours is closer to half that at $\sim 0.025$ dex.

The degree to which trends in the Galactic disk are azimuthally symmetric has the potential to provide interesting insight into the history of the disk. The stellar distribution across the Galaxy is not uniform, with in-situ non-axisymmetric features such as the Galactic bar and spiral arms containing higher stellar density than surrounding populations, particularly for young stars (e.g., Bland-Hawthorn & Gerhard, 2016; Reid et al., 2019; Khoperskov et al., 2020a). Additionally, chemical enrichment is strongly dependent on local conditions, and spiral density fluctuations can lead to measurable differences in a galaxy’s enrichment history across azimuth (Spitoni et al., 2019b).

To look more closely at the azimuthal symmetry of the disk, the right column of Figure 11 is a metallicity map of the disk identical to Figure 6 but displayed in polar coordinates. The spatial bins are sized $\Delta R=1$ kpc and $\Delta\theta=10^{\circ}$. As before, the stars are separated into rows based on their height above the Galactic plane. The corresponding panels to the left trace the median metallicity at each radius (y-axis) for different bins in azimuthal angle $\theta$ (point color), restricted to $130\leq\theta\leq 230$ deg, where there is reasonable coverage with radius. At each radius, the expected spread based on uncertainty of the median measurements is shown as the gray shaded region. The expected spread ($\sigma_{\rm[Fe/H]}$) is defined for a given radius as the sum of the uncertainty in the median for each individual $\theta$ bin ($\sigma_{\rm[Fe/H]}(R,\theta)$), divided by the number of bins with valid data ($N_{bins,\theta}$), whereas the uncertainty in the median for a given azimuth bin is the standard deviation in [Fe/H] divided by the square root of the number of stars in that bin:

Close to the Galactic plane ($|Z|\leq 0.5$ kpc, top row), the Galaxy has little metallicity variation across azimuth near the Solar neighborhood and outward ($R\geq 8$ kpc), although the region of azimuth covered by the observations decreases with increasing radius. Closer to the center of the Galaxy ($R\leq 5$) kpc, the median metallicity varies more with azimuthal angle $\theta$, with the spread possibly slightly exceeding the expected uncertainty. In the middle row higher above the Galactic plane ($0.5\leq|Z|\leq 1.0$ kpc), a similar trend is seen, where there is more spread with metallicity in azimuth near the center of the Galaxy. There does seem to be an asymmetry in the disk that follows the approximate location of the Galactic bar in these coordinates, with metal-rich stars preferentially residing in the Galactic bar. Higher above the Galactic plane ($1.0\leq|Z|\leq 2.0$ kpc, bottom row), the variation in azimuth can be attributed entirely to noise from low-number statistics, where the observed spread is all comparable or smaller than the expected spread.

Interactions with satellite galaxies and merger events can also perturb the disk in non-axisymmetric ways, such as warping the disk or introducing kinematic oscillations (e.g., Gómez et al., 2012; Cheng et al., 2020; Chrobáková et al., 2022). The restorative force from these perturbations is weaker in the outer Galaxy, so the presence of these features is generally expected to be more obvious at large radii. We find no significant azimuthal asymmetries in metallicity in the outer disk, but as radius increases, our sample covers less range in $\theta$. As such, we are unable to draw any conclusions on the Milky Way’s merger history from these metallicity maps alone.

Recent work has detected azimuthal variations in both the gas phase metallicity (Wenger et al., 2019) and stellar metallicity (Inno et al., 2019; Poggio et al., 2022; Hawkins, 2022) possibly corresponding to the Milky Way’s spiral arms. Measuring the radial metallicity gradient using Gaia data, Poggio et al. (2022) and Hawkins (2022) report azimuthal variation in the slope on the order of $0.02-0.1$ dex kpc^-1. If we measure the slope of the profiles in the left column of Figure 11 as a function of azimuthal angle $\theta$, our maximum difference between slopes for the vertical bin closest to the disk is 0.026, which is comparable to the lower end of the variation reported by Poggio et al. (2022), but may not be significant given the uncertainties in our data.

In Figure 12, we zoom into a region in the solar neighborhood, where the sample has higher number of stars making it possible to study the chemical distribution in more detail. The exact window used is outlined as the black rectangle in the top-right panel of Figure 11 for reference. Here, the spatial bins are sized $\Delta R=0.5$ kpc and $\Delta\theta=5^{\circ}$ (half the size of the bins in Figure 11), but calculated on a frequency of $\Delta R=0.1$ kpc and $\Delta\theta=1^{\circ}$ as a running median for smoothing. The approximate location of the nearby spiral arms from Reid et al. (2019) are plotted as colored lines. There does seem to be a bit of coherent structure signifying that the median metallicity is not symmetric in azimuth, but it does not obviously follow the spiral arms.

The radial and vertical metallicity gradients in the Milky Way disk have been well-documented observationally, with stars near the center of the Galaxy exhibiting higher metallicity than those at large radii and higher $Z$ (e.g., Hartkopf & Yoss, 1982; Cheng et al., 2012; Carrell et al., 2012; Anders et al., 2014; Schlesinger et al., 2014; Hayden et al., 2014; Frankel et al., 2019; Katz et al., 2021; Gaia Collaboration et al., 2022). Such a trend is predicted by “inside-out” disk formation models, where stars in the central regions of the Galaxy form earlier on in the Galaxy’s history, and the disk subsequently grows outward over time, the global star formation rate consistently decreasing with radius (e.g., Eggen et al., 1962; Larson, 1976; Matteucci & Francois, 1989; Kobayashi & Nakasato, 2011; Minchev et al., 2015). However, radial migration could complicate this interpretation because it flattens gradients over time as stars move away from their birth location (e.g., Sellwood & Binney, 2002; Roškar et al., 2008; Wang & Zhao, 2013; Hayden et al., 2015; Mackereth et al., 2017; Frankel et al., 2018, 2020; Vickers et al., 2021; Lian et al., 2022b). Without radial migration, gradients are predicted to steepen with lookback time, with the oldest stars having the steepest slope and being more centrally concentrated as a result of the inside-out growth of the Galaxy (e.g., Matteucci & Francois, 1989; Bird et al., 2013; Pilkington & Gibson, 2012; Gibson et al., 2013; Mollá et al., 2018). Radially dependent outflow efficiencies can also have strong impact on the radial gradients (Johnson et al., 2021), as can radial gas flows within the disk (e.g., Bilitewski & Schönrich, 2012).

We measure the median metallicity profile for the disk by first separating stars into bins of $Z$, and then calculating a running median [Fe/H] for each bin of $N=200$ data points with an overlap of 50%, sorted by Galactocentric radii. We do this for the total stellar population and repeat the analysis for the high-$\alpha$ and low-$\alpha$ samples described in Section 2.2 separately. The resulting radial median metallicity profiles are shown in Figure 13 and quantified in Table 1.

The total sample (Fig. 13 top row, Table 1) shows a negative metallicity gradient close to the Galactic plane, which flattens out at small radii. Moving above the plane, the slope of the gradient flattens until it becomes slightly positive at $|Z|>1.6$ kpc.

The low-$\alpha$ disk (Fig. 13 middle row) shows a steep metallicity gradient everywhere, notably missing the flattening in the inner Galaxy seen in the total population. The low-$\alpha$ disk’s metallicity gradient flattens with height $Z$ much like the total population.

The high-$\alpha$ disk (Fig. 13 bottom row) exhibits a much flatter or slightly positive metallicity profile whose slope does not change significantly with $Z$. The high-$\alpha$ sequence effectively ends at $R\gtrsim 10$ kpc (shown previously in Figure 8, meaning there are not enough high-$\alpha$ stars in the outer Galaxy to constrain the metallicity profile past $R\gtrsim 10$ kpc.

The total disk looks like the high-$\alpha$ profile near the center of the Galaxy, and matches the low-$\alpha$ profile in the outer Galaxy. This is due to the relative weights between these two populations at different locations: as shown in Section 3.2, the inner region of the Galaxy is dominated by the high-$\alpha$ sequence, and the outer region is mostly low-$\alpha$ stars.

For each median metallicity profile, we quantify the gradient by fitting a straight line to stars with Galactocentric radius $R\geq 7$ kpc, where the profile reasonably approximates a single line. The best-fit slope for each profile is shown in Figure 14 against height above the plane $Z$, and tabulated in Table 1. The total population and the low-$\alpha$ population have steep negative profiles in the outer Galaxy, which approach zero as it moves above the plane. The high-$\alpha$ slope is close to zero everywhere. Note that if we change the definition of our measured gradient and instead fit without the radial limit of ($R\geq 7$ kpc), the high-$\alpha$ population show a slight positive gradient, consistent with other studies (e.g., Vickers et al., 2021).

These results are generally consistent with previous results from a variety of methodologies. We measure a slope of $-0.056\pm 0.001$ dex kpc^-1 for the total population close to the Galactic plane ($|Z|\leq 0.25$ kpc). Using previous data releases of APOGEE, Feuillet et al. (2019) measured the slope of the low-$\alpha$ metallicity gradient to be $-0.059\pm 0.010$ dex kpc^-1. Using open clusters as tracers, Donor et al. (2020) measured a radial gradient of $-0.068\pm 0.001$ dex kpc^-1 with APOGEE DR16 data, and Myers et al. (2022) measured $-0.073\pm 0.002$ dex kpc^-1 with DR17. Using Gaia DR3 data, Gaia Collaboration et al. (2022) measured a slope of $-0.056\pm 0.007$ dex kpc^-1 for their bin closest to the Galactic plane. Additional studies use Cepheid stars as tracers and find similar results, with Genovali et al. (2014) reporting $-0.060\pm 0.002$ dex kpc^-1, and Lemasle et al. (2018) reporting $-0.045\pm 0.007$ dex kpc^-1. The large sample size, distance range, and high precision of the APOGEE DR17 sample enable us to map the radial, vertical, and $\alpha$-dependence of the metallicity gradient in unprecedented detail. Differences between tracer populations used may explain the slight differences among previous results — e.g., Cepheids tend to be young stars that are more concentrated to the mid-plane, which alone results in a steeper vertical gradient.

These findings are also generally consistent with Hayden et al. (2014), who measured the metallicity gradients as a function of position ($R$ and $Z$) in the disk and documented a relatively flat gradient near the center of the Galaxy, a steeper gradient farther out in the disk, and generally flat gradients for the high-$\alpha$ population (Hayden et al., 2014, Table 2). Our data set is significantly larger than that of Hayden et al. (2014), with the inclusion of Southern Hemisphere observations, which leads to better spatial coverage in both the inner and outer regions of the Galaxy and less sensitivity to potential systematics. This may be why our radial gradient $-0.056\pm 0.001$ dex kpc^-1 is slightly shallower than the $-0.073\pm 0.003$ dex kpc^-1 from Hayden et al. (2014) for a comparable spatial zone.

Numerical simulations from Rahimi et al. (2014) reproduce similar gradient trends where the radial metallicity gradients flatten with increasing $Z$. Notably, they attribute the slight positive gradient of the Milky Way’s thick disk to the flaring of younger populations at large radii. The vertical flaring of the mono-age, mono-abundance disk has been documented extensively in other studies as well (e.g., Minchev et al., 2015; Bovy et al., 2016; Mackereth et al., 2017; Lian et al., 2022a; Gaia Collaboration et al., 2022; Robin et al., 2022)).

The models of Johnson et al. (2021), which incorporate radial and vertical redistribution of stars based on a hydrodynamical cosmological simulation of a Milky Way-like galaxy, also show a radial metallicity gradient that flattens with increasing $|Z|$.

The vertical median metallicity profile of the disk is shown in Figure 15 and Table 2, calculated in the same way as the radial gradients. The best slopes for all stars $|Z|<2$ kpc are shown in Figure 16 and Table 2. Consistent with Gaia Collaboration et al. (2022), we find the gradients in the disk to be vertically symmetric. The slope measured above the disk ($Z>0$) does not differ significantly from the slope measured below the disk ($Z<0$). The total population has a steep negative vertical gradient close to the center of the Galaxy ($R<5$ kpc), which flattens moving out in radius. This is generally true for the low-$\alpha$ population as well, although the innermost parts of the Galaxy ($R<2$ kpc), show a very flat profile, possibly due to the bulge or bar. The high-$\alpha$ population has a shallow negative gradient everywhere, which does not significantly change with radius. Beyond $R\gtrsim 10$ kpc, the vertical gradient for the total, high-$\alpha$ and low-$\alpha$ populations are all close to zero.

This is also consistent with previous studies (e.g., Hayden et al., 2014, Table 1), where the vertical metallicity gradient approaches 0 as radius increases. Near the solar neighborhood, we report a vertical median metallicity gradient of $-0.315\pm 0.009$ dex kpc^-1. Hayden et al. (2014) report a gradient of $-0.31\pm 0.01$ dex kpc^-1 in a comparable spatial zone. For the thick disk, Carrell et al. (2012) measured a vertical metallicity gradient of $-0.113\pm 0.010$, using stars close to the solar neighborhood with heights $1\leq|Z|\leq 3$, consistent with our measurement of $-0.112\pm 0.014$ near the solar neighborhood for the high-$\alpha$ population.

Examining the Galaxy’s metallicity profile as a function of age provides a direct link to the evolution history of the disk. Figure 17 shows the radial (top panel) and vertical (bottom panel) metallicity profile for the low-$\alpha$ disk only, separated into bins of stellar age. In both cases, the slope of the profile flattens with increasing age. This trend is also seen in Anders et al. (2023), using a similar APOGEE sample with different age and distance estimates. The flattening of the metallicity gradients with age is the opposite of what is predicted by a pure inside-out growth of the Galaxy (e.g., Matteucci & Francois, 1989; Bird et al., 2013), where the gradient in the interstellar medium is expected to flatten out over time (e.g., Pilkington & Gibson, 2012; Gibson et al., 2013; Mollá et al., 2018). The opposite trend, seen here, is commonly attributed to be a signature of radial migration (e.g., Wang & Zhao, 2013; Magrini et al., 2016; Minchev et al., 2018; Vickers et al., 2021; Anders et al., 2023). However, an alternate explanation is also presented in Chiappini et al. (2001), where a disk formed from pre-enriched gas starts with an initially flat metallicity gradient that steepens over time.

While the radial and vertical median metallicity gradients in the disk reveal interesting general trends, more insights can be gleaned from the full metallicity distribution function (MDF) at different locations in the disk. Specifically, the spread and shape of the underlying MDF can be crucial for characterizing the complex history of the disk more accurately.

Figure 18 demonstrates how the metallicity distribution function varies with radius for samples at different vertical slices, from closest to the Galactic plane (top panel) to furthest beyond (bottom panel). Every third row is annotated with tick marks denoting the 25th, 50th (median), and 75th percentile of that row’s distribution, for the total sample (white), low-$\alpha$ (blue), and high-$\alpha$ (red) samples separately. The diamond point is the peak (or mode) of the distribution. At all heights, the characteristic metallicity (whether median or mode) decreases with radius for the total and low-$\alpha$ populations, but stays roughly constant for the high-$\alpha$ population.

At smaller radii, there is little overlap between the low-$\alpha$ and high-$\alpha$ MDF. In fact, the high-$\alpha$ population alone is what creates the metal-poor tail of the total MDF. Just outside the solar neighborhood ($8<R<12$ kpc), the MDF of the high-$\alpha$ and low-$\alpha$ populations overlap chemically. Moving above the plane, this overlap starts closer to the center of the Galaxy, around $R=5$ kpc at $|Z|>1.0$ kpc.

The shape of the MDF changes as a function of radius; close to the center of the Galaxy, the total and low-$\alpha$ distribution is heavily skewed towards lower metallicities (the peak trends right of the median), and in the outer Galaxy the distribution is skewed towards higher metallicities (the peak is left of the median). This is consistent with the trend seen in previous studies (e.g., Anders et al., 2014; Kordopatis et al., 2015b; Hayden et al., 2015; Loebman et al., 2016; Katz et al., 2021).

The high-$\alpha$ MDF is consistently broader than the low-$\alpha$, has a shallower characteristic gradient, and shows less of a skewness trend with radius; the peaks (diamond points) are closer to the median (center tick mark) in general. Stars in the low-$\alpha$ sample transition from being negatively skewed in the inner Galaxy to positively skewed in the outer Galaxy (as shown in Figure 19). However, this trend is not seen as strongly in the high-$\alpha$ disk, even at similar Galactic heights as the thin disk. Because the high-$\alpha$ population is generally older, this could imply that the high-$\alpha$ sample is more “mixed” vertically, meaning the birth location of stars tend to be farther away from their present-day location, largely because they have had more time to migrate radially.

Kordopatis et al. (2015b) also observed an overabundance of metal-rich stars in the solar neighborhood using data from the RAVE survey, and link it to radial migration. Hayden et al. (2015) observed this same trend in APOGEE and present a simple model to show that radial migration could explain the change of the MDF shape with radius, because more stars migrate outward from the inner disk than vice versa. Loebman et al. (2016) and Johnson et al. (2021) show that this explanation succeeds quantitatively in models with realistic radial migration from cosmological simulations. These studies also show that the excess metal-rich tail in the MDF disappears at high $|Z|$, in agreement with the results of Figure 18 and Hayden et al. (2015).

Some simple statistics can be measured to more fully characterize the MDF as a function of radius, as shown in Figure 19 for the case of the low-$\alpha$ disk close to the Galactic plane (equivalent to the top panel of Figure 18). The top right plot shows three definitions for a characteristic value for [Fe/H]; the mean, median, and peak of each distribution as a function of radius. While these values are similar, they are not identical, meaning the measured metallicity gradient of the disk will depend on the parameterization of [Fe/H] chosen. In the inner disk, the peak [Fe/H] is up to 0.1 dex higher in metallicity than the mean and median, due to the distributions being skewed metal-rich. In the outer disk, the peak is preferentially more metal-poor. Therefore, a metallicity gradient measured using the peak metallicity as a tracer will have a steeper slope than a gradient measured with the mean metallicity for the same group of stars.

The middle right panel in Figure 19 quantifies the spread of each distribution with radius, as total standard deviation $\sigma$. The outer regions of the disk ($R>10$ kpc) are characterized by narrower distributions with less overall spread, whereas the inner disk MDFs span a larger range of metallicities.

The bottom right panel quantifies the skewness of each distribution with radius. The MDF in the inner regions of the disk it is negatively skewed, and in the outer regions it is positively skewed, quantifying the trend seen earlier Figure 18.

The present-day distribution of stellar ages can act as an interesting snapshot as to what the Milky Way might have looked like at different points in time, while also documenting how stars might move and migrate away from their radius of birth.

The radial median age profile of the Milky Way is presented in Figure 20, identical to the way the metallicity gradients in Section 3.4 were calculated. The best-fit slope for each profile is shown in Figure 21, once again calculated only using stars with $R\geq 7$ kpc where the profile reasonably approximates a straight line.

The median age of the high-$\alpha$ population is generally the same everywhere, with a slope close to 0 at any height above the plane. We note that the median age observed here, $\sim 8.5$ Gyr, is likely too young due to selection effects and the limitations of our age sample, including the metallicity cut of [Fe/H] $\geq-0.7$ dex discussed in Section 2.4.

The total and low-$\alpha$ stellar populations have a negative radial median age gradient in the outer Galaxy, while in the inner Galaxy the profile flattens out. The measured slope is flattest close to the disk ($z=0$), and becomes negative moving above the plane. The total and low-$\alpha$ stellar populations have a vertical gradient that varies with radius as well. The vertical gradient is steeper closer to the center of the Galaxy, and flatter at large radii, similar to the vertical metallicity gradient. Unlike the high-$\alpha$ sample, the majority of stars in the low-$\alpha$ disk have metallicity [Fe/H] $\geq-0.7$ and are therefore more robust against potential biases caused by the sample selection (see Figure A.1).

As before with metallicities, more information lies in the shape of the age distribution function at different locations in the Galaxy rather than the gradient alone. Figure 22 depicts the age distribution function (ADF) as a function of Galactocentric radius for different heights in the disk. Close to the Galactic plane (left panel), the peak age of the total population gradually declines, peaking around 7.5 Gyr near the Galactic center and 3 Gyr near the solar radius and outward. The spread of the ADF is fairly broad, spanning up to 5 Gyr at all radii. In the outer Galaxy, the ADF is preferentially skewed toward older ages. Katz et al. (2021) found the ADF skewed towards younger ages in the inner Galaxy, and skewed towards older ages in the outer Galaxy, but they limited their investigation to just the low-$\alpha$, thin disk sample.

Farther above the Galactic plane (middle and bottom panels of Fig. 22), the profile does not show a single gradient, but rather has a slight positive age gradient until $R\sim 5$ kpc which transitions into a negative gradient at larger $R$. After $R\sim 12$ kpc, the gradient flattens out. We caution that this flattening of the gradient at large $R$ may be artificially induced by the lack of [Fe/H]$<-0.7$ stars in our sample when dealing with ages imposed by the sample cuts described in Section 2.4.

Separating further into the low-$\alpha$ and high-$\alpha$ samples reveals slightly different trends. Close to the Galactic plane (top panel), and in the inner Galaxy, there is some minimal overlap between the low-$\alpha$ and high-$\alpha$ samples, but near the solar neighborhood and outward, the ADF is more bimodal and there is more separation between the high-$\alpha$ and low-$\alpha$ ADF. At all radii, the low-$\alpha$ sample has a narrower distribution, and transitions from being skewed toward younger ages in the inner Galaxy to being skewed toward older ages in the outer Galaxy. This is similar behavior to the MDF in Figure 18, and also seen by Katz et al. (2021). Farther above the plane, there is more overlap between the low-$\alpha$ and high-$\alpha$ samples in the inner Galaxy, but there is still little overlap at larger radii.

Despite the similarities in these representations of the ADF to the corresponding ones for MDFs, we caution that the age uncertainties (typically $\sim 0.15$ dex) are significant relative to the total spread, which is not the case for the [Fe/H] measurements.

The relation between stellar age and metallicity has long been sought after to help constrain chemical evolution models (e.g., Twarog, 1980; Edvardsson et al., 1993). In a simple “closed box” system, the metallicity of stars should increase over time, as each generation of stars enriches the interstellar gas from which subsequent generations are born. The actual scenario is much more complex, depending on gas inflow and outflow rates, supernovae yields, stellar migration, and the positionally variable star formation history of the Galaxy. Observations of the age-metallicity relation in the Milky Way include the effects of all of these processes, and they therefore provide a powerful constraint for chemical evolution models to reproduce across different locations in the Galaxy.

Previous studies have found significant scatter in the age-metallicity relation (AMR) near the solar neighborhood, where stars with a single age span a wide range of metallicities, which cannot be attributed to observational errors alone (e.g., Casagrande et al., 2011; Bergemann et al., 2014; Aguirre et al., 2018; Lin et al., 2018; Grieves et al., 2018; Sahlholdt et al., 2022; Xiang & Rix, 2022; Anders et al., 2023). The AMR also varies with Galactic location, making it difficult to constrain a single relation that fits the whole disk (e.g., Hasselquist et al., 2019; Feuillet et al., 2019; Casamiquela et al., 2021; Lian et al., 2022b; Anders et al., 2023).

The AMR for our sample of stars in the Milky Way disk is shown in Figure 23. The data are split into bins of radius (columns), and height (rows), to demonstrate how the age-metallicity relation varies across Galactic location. Points are colored by $\alpha$-element abundances, which are known to be correlated with age, although the exact trend depends on Galactic position (e.g., Haywood et al., 2013; Nissen, 2015; Bedell et al., 2018; Feuillet et al., 2018; Hasselquist et al., 2019; Lian et al., 2022b, see also Section 3.9). The median trend across different bins in age is tracked as black square points to help guide the eye. Note that the x-axis (stellar age) has been reversed so that older stars are on the left, better expressing the forward flow of time. The x-axis is also presented in log space, which is more representative of the uncertainties in our age estimates (Section 2.4). One can also read each panel rotated 90^∘ to see the distribution of log(age) at fixed [Fe/H].

Consistent with previous studies, there is significant scatter around the age-metallicity relation near the solar neighborhood. There is a slight gradient in $\alpha$-element abundances within that spread; nearly everywhere in the Galaxy, higher-metallicity stars are relatively more $\alpha$-poor.

The age-metallicity relation has less scatter moving towards the inner Galaxy (left column), and above the plane (top row). If the spread in this relation is due to the radial migration of stars, this implies that the in-situ age-metallicity relation of the inner galaxy has been more preserved, and less contaminated by migrated stars; or in other words, more stars migrate outwards than inwards in the Galaxy. This is perhaps not surprising; dynamically, stars in the inner Galaxy migrate outwards, and stars in the outer Galaxy migrate inwards (e.g., Sellwood & Binney, 2002). The inner Galaxy is denser than the outer disk, so if one assumes the same rate of migration across all radii, more stars would migrate outwards simply because more stars start in the inner Galaxy. The age-metallicity relation shown in Sahlholdt et al. (2022) also shows a lower dispersion for the inner regions of the disk (their “Pop C” sample) using GALAH data.

The AMR is steepest in the inner Galaxy, and flattens moving out in radius. In all panels, the age-metallicity relation flattens out at young stellar ages, perhaps indicating that chemical equilibrium has been reached, where the inflowing gas dilutes the ISM at the same rate as it is being enriched (e.g., Dalcanton, 2007; Finlator & Davé, 2008; Weinberg et al., 2017). In the inner disk, the equilibrium metallicity is higher than the equilibrium metallicity reached in the outer disk and moving away from the midplane. Equilibrium seems to have been reached sooner (at older stellar age) in the outer disk than the inner disk.

There is a notable inversion of the AMR at large radii in the Galaxy, where older stars trend more metal-rich than the younger stars in the sample. This has been reported before in previous studies (e.g., Anders et al., 2014; Feuillet et al., 2018; Hasselquist et al., 2019; Lian et al., 2022b). The cosmological simulations of Lu et al. (2022) explore the possible origins of such an inversion, suggesting that it could be the signature of interactions with a satellite galaxy like the Sagittarius dwarf, and radial migration notably widens the apex.

A recent study from Xiang & Rix (2022) documented a disjointed age-metallicity relation for the sum of the total disk, and suggested that the two-infall scenario or a major merger is responsible. In our results, we see some evidence of bimodality in Figure 24, which is a contour map of the point density distribution in Figure 23. In the inner Galaxy ($R<3$ kpc), the age-metallicity relation has two separate peaks.

While no clear correlation between age and metallicity relation exists near the solar neighborhood, better correlation has historically been found between age and $\alpha$-element abundances (e.g., da Silva et al., 2012; Haywood et al., 2013; Bensby et al., 2014), leading some to use $\alpha$-elements as a chemical clock, substituting their abundances when stellar ages are not readily available. Even so, the [Mg/Fe]-age correlation has been found to vary across the Milky Way’s disk (e.g., Aguirre et al., 2018; Feuillet et al., 2018; Katz et al., 2021; Vázquez et al., 2022), extending the metaphor to imply that chemical clocks run in chemical “time zones” throughout the Galaxy. Reproducing this variation has been considered a strong constraint on chemical evolution models (e.g., Haywood et al., 2013; Spitoni et al., 2019a; Johnson et al., 2021).

Figure 25 shows the relation between age and [Mg/Fe] as a contour plot at various locations throughout the disk. The distribution is double-peaked in nearly all panels, with the older, low-$\alpha$ population most prevalent above the disk and in the inner Galaxy. Near the solar neighborhood and outward, the relation of the low-$\alpha$ sequence is relatively flat, with a large spread in ages corresponding to a small range of [Mg/Fe] abundances. This is consistent with the findings of Haywood et al. (2013) and Feuillet et al. (2018). In the inner Galaxy, a small range in [Mg/Fe] abundances corresponds more tightly with a smaller range in age, a phenomenon that applies, though differently, to each the high-$\alpha$ and low-$\alpha$ groups. As in Haywood et al. (2013), we observe some age overlap between the two sequences, implying that the low-$\alpha$ sequence in the outer disk began forming stars while the high-$\alpha$ disk was concurrently still forming stars in the center of the Galaxy.

For the low-$\alpha$ sequence, there is some evolution with Galactic position. The low-$\alpha$ stars are older and more $\alpha$-poor near the center of the Galaxy. In the outer Galaxy, the sequence is more $\alpha$-enhanced and generally younger, although covering a larger spread in ages. Above the plane ($Z>1$ kpc), the low-$\alpha$ sequence is more $\alpha$-enhanced.

The high-$\alpha$ sequence stays generally in the same location on this diagram regardless of position in the Galaxy. This is similar to the trend seen in Figure 8 and 10, where the locus of the low-$\alpha$ sequence changes significantly with Galactic position while the high-$\alpha$ sequence stays largely in the same location.

In the $R=9-12$ kpc and low-$|Z|$ zones, the log(age) distribution is bimodal even within the low-$\alpha$ sequence. This could be evidence for a three-phase star formation history. Sahlholdt et al. (2022) report a similar distribution in the age-metallicity relation of their “Pop A” sample, which probes a similar location in the disk. The younger peak in log(age) appears similar to the recent starburst 2-3 Gyr ago detected independently in Isern (2019) and Mor et al. (2019), although we note that this is the first time to our knowledge that it has been detected in APOGEE data. This recent starburst is thought to be linked with the Sagittarius dwarf spheroidal galaxy’s most recent perigalactic passage through the Milky Way’s disk (Ruiz-Lara et al., 2020; Laporte et al., 2019; Antoja et al., 2020). The uncertainties in our age estimates are not negligible, but are likely not responsible for this bimodality. Larger uncertainties would blur out the distribution and decrease the observed bimodality.

We remind the reader that using different estimates for stellar ages can change our results. Alternate versions of this figure using different age catalogs are shown in Appendix A to emphasize this point.

If a group of stars was born together in the same location and at the same time, they should have identical chemical abundances ([Fe/H] and [Mg/Fe] in our case). This is the basic assumption behind “chemical tagging”, used to identify stellar siblings that have been redistributed throughout the Galaxy despite being born together (e.g., Freeman & Bland-Hawthorn, 2002; Hawkins et al., 2015; Ting et al., 2015; Price-Jones et al., 2020; Buder et al., 2021). Under this assumption, we can potentially track radial migration and the spatial redistribution of a stellar population over time, as well as look at the enrichment history of an area of the Galaxy for fixed metallicity.

Figure 26 shows this evolution for stars close to solar metallicity ($-0.1\leq$ [Fe/H] $\leq 0.1$ dex). The present-day radial distribution (x-axis) and [Mg/Fe] abundance (y-axis) is shown for different bins in stellar age (line color) for both the low-$\alpha$ (solid lines) and high-$\alpha$ populations as the contour containing 75% of all points in that bin. Stars at the same age and metallicity should all have same [Mg/Fe] abundance. If stars did not move from where they were born, we would expect to see a tight clump in this space. If stars migrate significantly over time, the shape of the clump should spread out to a broader range of radii, but keep the same [Mg/Fe] abundances (or a “flat slope” in [Mg/Fe]). If the slopes were not flat, it may be indicative of different enrichment histories for different parts of the Galaxy; suggesting a violation of the assumption that stars with the same metallicity and age were born roughly in the same place.

The young, $\alpha$-poor stars in Figure 26 currently reside in a relatively confined clump in radius and [Mg/Fe] as expected ($6\lesssim R\lesssim 11$ kpc). As stellar age increases, the shape of the clump widens to cover a broader range in radius ($3\lesssim R\lesssim 12$ kpc) while the [Mg/Fe] abundance remains confined. The high-$\alpha$ sequence does not show this evolution in radial width with time, but notably does not include enough young stars to properly trace this. As shown in Figure 9, high-$\alpha$ stars tend to be old.

As stellar age decreases, the median [Mg/Fe] value of each subpopulation decreases for both the low-$\alpha$ (solid lines) and high-$\alpha$ (dashed lines) populations. In the low-$\alpha$ population, the median value of [Mg/Fe] decreases from [Mg/Fe] = 0.05 in the oldest age bin ($10\leq\log(\textrm{age})\leq 10.1$) to [Mg/Fe]$=0.00$ in the youngest age bin ($9\leq\log(\textrm{age})\leq 9.1$). This tracks the chemical evolution of a location in the Galaxy, as Type Ia supernovae begin “diluting” the interstellar medium with iron, thereby decreasing the overall [Mg/Fe] ratio over time.

The “superposition” class of evolution models (e.g., Schönrich & Binney, 2009a, b; Minchev et al., 2013, 2014, 2017; Johnson et al., 2021) explain the observed chemical bimodality in the solar vicinity as the superposition of evolutionary tracks for stars born at different Galactocentric radii, with the stars having reached their present-day location in the solar neighborhood through radial migration. Several predictions made by these superposition chemical evolution models are seen in our results.

The metallicity gradient flattening with age (Figure 17) is predicted by radial migration (e.g., Sellwood & Binney, 2002; Roškar et al., 2008; Wang & Zhao, 2013; Hayden et al., 2015; Mackereth et al., 2017; Frankel et al., 2018, 2020; Vickers et al., 2021; Lian et al., 2022b). If stars are formed in-situ with a steep metallicity gradient, that gradient will flatten over time as metal-rich inner Galaxy stars migrate outwards and metal-poor outer Galaxy stars migrate inwards, skewing the metallicity distribution at either end of the Galaxy. The older stars in our sample show a flatter gradient than the younger stars; in agreement with this scenario.

The shape of the metallicity distribution function at different locations in the Galaxy (Figures 18 and 19), specifically the skewness or asymmetry of the MDF, can be a sign of radial migration if a population of stars has a metal-rich tail (e.g., Hayden et al., 2015; Kordopatis et al., 2015b; Loebman et al., 2016)). Our data show a skewed MDF in Figure 19, where the inner region of the disk is negatively skewed, and the outer region of the disk is positively skewed. This trend is commonly attributed to radial migration, whereby migrating stars become the metal-rich tails in the MDF at different locations (e.g., Schönrich & Binney, 2009a, b; Roškar et al., 2008; Kordopatis et al., 2015b; Loebman et al., 2016; Johnson et al., 2021). The metal-poor tail in the inner Galaxy’s MDF is likely attributed to a spread in ages between the stars, as is predicted by even closed-box chemical evolution tracks when any given location in the Galaxy becomes enriched over time (e.g., Romano & Starkenburg, 2013; Vincenzo et al., 2014; Weinberg et al., 2017; Toyouchi & Chiba, 2018). The metal-rich tail of the outer Galaxy’s MDF is more difficult to explain with a traditional chemical enrichment track, which leaves radial migration as the most likely culprit.

The inversion in skewness in the MDF occurs around $R=9.4$ kpc in our data. This could be linked to the Outer Lindblad Resonance (OLR) of the Milky Way disk; a resonance with the Galactic bar driving different dynamical effects throughout the disk, one of which is radial migration (Halle et al., 2015; Michtchenko et al., 2016; Dias et al., 2019; Khoperskov et al., 2020b). Using Gaia data, Khoperskov et al. (2020a) estimates the OLR to be located at a Galactocentric radius of around 9 kpc. Khoperskov et al. (2020b) uses a high-resolution N-body simulation to investigate the relationship between the OLR and radial migration, and find that stars from the inner Galaxy migrating outwards become “trapped” in the OLR. When the rotation period of the bar slows down, those stars can escape and migrate farther out. The trapping effect of the OLR could also explain the build-up of metal rich stars at $R\sim 9$ kpc in Figures 5 and 6.

We also see signs of radial migration in the age-metallicity relation shown in Figure 23. Around the solar neighborhood, there is significant scatter about the age-metallicity trend. If the interstellar medium is always well-mixed, the large spread in metallicity for stars of a given age must mean that some of these stars were not born at their present-day location. This is another consequence of radial migration predicted by the superposition class of models, which explain the spread by emphasizing the difference between the present-day locations of stars and their birth radii (e.g., Schönrich & Binney, 2009a; Wang & Zhao, 2013; Minchev et al., 2013; Lian et al., 2022b). While any location in the Galaxy should start with a tight age-metallicity relation, migration will blur the present-day relation as metal-rich stars from the inner disk move outward and contaminate the more metal-poor outer disk. Radial migration is most efficient in the plane of the disk, therefore these models predict less spread in the age-metallicity relation at larger vertical heights, an effect also seen in our data.

The turnover in the age-metallicity relation, seen in Figure 23 and in Hasselquist et al. (2019), can also be explained through radial migration. The older, metal-rich stars were likely formed in the inner Galaxy, and migrated outwards to where they are found today, contaminating the age-metallicity relation. An alternate explanation could be the dilution of the ISM from pristine gas infall, lowering the metallicity of a previously-enriched area of the Galaxy (e.g., Spitoni et al., 2019a; Lu et al., 2022). However, in this scenario, it is predicted that both the high-$\alpha$ and low-$\alpha$ tracks in metallicity would still decrease with stellar age, with the post-dilution low-$\alpha$ track beginning at a lower metallicity than the high-$\alpha$ track at the same time (e.g., Spitoni et al., 2019a). That predicted trend is not obvious in our data, and the large spread in metallicity at a given age favors a radial migration explanation. Lu et al. (2022) further explore the origin of such a turnover using cosmological simulations, and find that even when the turnover in the age-metallicity relation can be directly linked to the infall of a satellite galaxy, radial migration can widen the shape of the peak.

Vertical motions should also be considered. In the age-$\alpha$ relation of Figure 25, a population of stars with the same age and in the same present-day location can have a large spread in [Mg/Fe], most dramatically seen in the inner Galaxy ($R<3$ kpc). This violates the underlying assumption mentioned earlier, meaning stars must have moved around to create that spread. However, in this case, radial migration is an unsatisfactory explanation, as stars from the outer Galaxy migrating inward is expected to be a less frequent occurrence than the other way around, simply a consequence from the density profile of the Galaxy (e.g., Sellwood & Binney, 2002). Instead, the vertical motion of stars could result in the observed spread if stars born above the mid plane are currently found near $Z=0$. Even stars formed above the plane will inevitably have vertical motions that cause their orbits to cross the plane, meaning this could be a natural consequence of a star’s vertical orbit. Related vertical motions could be linked with dynamical heating (e.g., Spitzer & Schwarzschild, 1951; Barbanis & Woltjer, 1967; Lacey, 1984; Mackereth et al., 2019) or the “upside-down” formation of the disk (e.g., Toth & Ostriker, 1992; Quinn et al., 1993; Hänninen & Flynn, 2002; Brook et al., 2004; Freudenburg et al., 2017; Bird et al., 2013, 2021). Using a quantity like guiding radius ($R_{\rm guide}$) and maximum height $|Z_{\rm max}|$ calculated from parameterized stellar orbits, instead of the present-day $R$ and $|Z|$ we use here, may remove contamination by thick disk stars currently “passing through” the thin disk from these figures. Some studies, including Boeche et al. (2013), Katz et al. (2021), and Spitoni et al. (2022a), have looked at these quantities, and found similar overall trends.

The evolution of radial distribution with stellar age seen in Figure 26 is yet more evidence for radial migration. The youngest, low-$\alpha$ stars indicate that for a population of fixed age, [Fe/H], and [Mg/Fe], they are expected to be born at a similar radius in the Galaxy. As stellar age increases, the population redistributes into a larger range of $R$, showing that stars migrate radially over time. The high-$\alpha$ stars, which cover a similar range in Galactic radius despite a wide window in age ($5\leq$ age $\leq 13$ Gyr), may suggest an “upper limit” on the efficiency of radial migration and the timescales over which stars can migrate on average; see e.g., Frankel et al. (2020) and Lian et al. (2022b) for such an analysis. Similarly, the relative distribution of stars within the low-$\alpha$ contours may hint at the importance of direction in radial migration (i.e., what fraction of stars migrate outwards instead of inwards), although such a discussion is beyond the scope of this work.

This wide variety of results suggests that significant stellar migration occurs in the Milky Way disk, most influencing the trends seen close to the Galactic plane and at larger radii. However, there are open questions remaining about the nature of the inner disk, where an apparent bimodality exists that is not easily explained by migration models.

The “two-infall” class of evolution models (e.g., Chiappini et al., 1997, 2001; Lian et al., 2020a, b, c; Spitoni et al., 2019a, 2020, 2021, 2022a) suggest that the two sequences in [Mg/Fe]-[Fe/H] space formed sequentially in time. Under this scenario, the thick disk was formed first during the initial collapse of the Galaxy, and after some time delay a second infall of gas fed the creation of the thin disk. These models predict several of the observed trends in our results, most notably in areas where radial migration is not as efficient, including the inner Galaxy and at greater heights above the Galactic plane.

Figure 8 shows that the $\alpha$-bimodality persists throughout the majority of the disk. This is significant because while the “superposition” class of chemical evolution models can produce the low-$\alpha$ “sequence” and broad distribution of [Mg/Fe] in the solar neighborhood using only radial migration (e.g., Schönrich & Binney, 2009a, b; Minchev et al., 2013, 2014; Nidever et al., 2014; Sharma et al., 2021a; Johnson et al., 2021), radial migration is known to be most efficient close to the Galactic plane and in the outer disk. Therefore, the bimodality in other parts of the Galaxy is more difficult to explain with radial migration alone. Previous studies by Freudenburg et al. (2017) and Zasowski et al. (2019) report that the shape of the MDF and the [Mg/Fe]-[Fe/H] trends seen in the inner disk ($3<R<5$ kpc) could be modeled well using a single evolutionary track in an “upside-down” disk formation model.

Our results are consistent with this in the [Mg/Fe]-[Fe/H] realm, but when expanded to ages, we show that two distinct tracks of [Mg/Fe] are observed even in the inner Galaxy. This is also reported in Queiroz et al. (2021).

Near the solar neighborhood, radial migration models predict bimodality by explaining the high-$\alpha$ sequence as contaminants from the inner Galaxy; due to the intrinsic density profile of the Galaxy, more stars are expected to migrate outwards than inwards, so the inner Galaxy should display less contamination from the low-$\alpha$ sequence. This is in contradiction with our findings, where the $\alpha$-bimodality persists throughout the majority of the disk. Johnson et al. (2021) also report that the bimodality reproduced by their superposition model is weaker than the observed bimodality in the Milky Way; in particular, the model overproduces intermediate-$\alpha$ stars compared to observations. Chen et al. (2022) find greater success in producing bimodality with radial migration. However, a two-infall (or multi-infall) model may be needed to explain the bimodality in the inner Galaxy.

The age-metallicity relation in Figure 23 produces some trends that are better explained by the two-infall model than by radial migration. Minchev et al. (2013, 2014) report that while the scatter around the age-metallicity relation can be attributed to radial migration, the overall slope is only weakly affected. In our data, consistent with other recent studies (Feuillet et al., 2019; Hasselquist et al., 2019), the slope of the age-metallicity relation varies significantly with Galactic radius. This is reproduced with a two-infall model, where the low-$\alpha$, post-infall disk has a shallower slope in the age-metallicity relation compared to the high-$\alpha$ population, due to the continuous inflow of gas diluting the disk that was not present during the formation of the original high-$\alpha$ disk (Spitoni et al., 2019a, 2020, 2021).

The apparent disk bimodality is not only observed in [Mg/Fe]-[Fe/H] chemistry, but also in the age-metallicity relation. Recent work by Xiang & Rix (2022) show a disjointed age-metallicity relation for the sum of the total disk, which is only possible to reproduce in a two-infall scenario or with a major merger event. In our study, the age-metallicity relation as a function of Galactic position (Fig. 23) also shows a possible bimodality, most apparent in the $0<R<3$ kpc range of the inner disk.

The age-[Mg/Fe] relation in Figure 25 is perhaps more convincing evidence for the two-infall model, with clear bimodality in the relation persisting across nearly the entire disk. In the outer disk ($9<R<15$), there appears to be an additional bimodality within the low-$\alpha$ sequence, suggesting a three-phase star history similar to that detected in Sahlholdt et al. (2022) using a sample of stars from the GALAH survey. The uncertainties in our age estimates are not negligible, but are likely not responsible for this bimodality. Larger data uncertainties would blur out the distribution and decrease the observed bimodality. Due to the age uncertainties, the “true” age bimodality in the Milky Way may be stronger than what is shown in our analysis. The peak of our “third” starburst is around 2-3 Gyr, consistent with the recent burst Isern (2019) and Mor et al. (2019) report, possibly linked to the most recent interaction with the Sagittarius dwarf spheroidal galaxy as it passed through the disk (e.g., Ruiz-Lara et al., 2020; Laporte et al., 2019; Antoja et al., 2020)

In the age-[Mg/Fe] relation, the age overlap between the high-$\alpha$ and low-$\alpha$ sequences are impossible to explain with a single evolutionary track. In this interpretation, we do caution that if the transition from high-$\alpha$ to low-$\alpha$ is fast, uncertainties in age determination could produce an artificial impression of age overlap. Haywood et al. (2013) also detects an overlap in age, and claim that the dichotomy in the solar neighborhood can be reproduced by the two-infall scenario, and that little to no radial migration is needed.

For the inner disk, the spread in [Mg/Fe] for a given age is likely due to the vertical motions of stars as discussed previously in Section 4.1. However, if the spread was purely from vertical blurring, a continuous spread would be expected, and not the bimodal distribution seen here.

As discussed in Appendix A, the choice of stellar age catalog does significantly change the appearance of the age-[Mg/Fe] relation across the disk, suggesting that stellar age estimates are not yet robust enough to draw strong conclusions from this relation. Nevertheless, several of the key features observed here (the age overlap and the bimodality within the low-$\alpha$ sequence) do persist across different age catalogs.

In summary, the bimodality in [Mg/Fe] and stellar ages persisting across the inner disk is not easily explained through radial migration, which is most efficient at larger radii. A multi-phase star formation history, such as those presented in the two-infall model, better predicts the trends observed in the inner Galaxy.

The “clumpy star formation” models (e.g., Clarke et al., 2019; Beraldo e Silva et al., 2020; Amarante et al., 2020) predict that the two sequences in $\alpha$-element abundances formed simultaneously but in different modes: the high-$\alpha$ sequence formed in rapidly-enriched gaseous clumps and the low-$\alpha$ formed in a less efficient smooth disk. This clumpy phase of early disk formation is predicted by numerical simulations (e.g., Bournaud et al., 2007) and often seen in observations of high redshift galaxies (e.g., Elmegreen et al., 2005).

One major result that the clump star formation models predict is the temporal overlap between the two $\alpha$ sequences: The low-$\alpha$ disk starts forming at the same time as the high-$\alpha$ sequence, meaning there should be some overlap in stellar ages between the two sequences. We see this overlap in Figure 23 and 25, where stars at around $\log$(age) $\sim 9.7$ Gyr span a significant range in both metallicity and $\alpha$-element abundances. This age overlap has also been observed in previous studies (e.g., Haywood et al., 2013; Hayden et al., 2017; Aguirre et al., 2018; Gent et al., 2022). Unfortunately, high uncertainties in age estimates could artificially produce this overlap.

The age distribution in the [Mg/Fe]-[Fe/H] plane shown in Figure 9 also possibly points to this formation scenario. As discussed in Section 3.2, the distribution of stellar ages within the low-$\alpha$ sequence suggests that it is not a single evolutionary sequence. Stellar age is more closely correlated with [Mg/Fe] instead of [Fe/H]. A strikingly similar trend is notably predicted under the clumpy formation scenario (Clarke et al., 2019, Fig. 9).

One potential avenue for further investigating the difference between the two-infall class of models and the clumpy star formation models lies not within the Milky Way but in other galaxies. If the two-infall model is true, a chemical bimodality would only be present in galaxies that experienced significant gas infall both at early and late times, meaning it would be a rare phenomenon only affecting approximately 5% of galaxies with comparable mass to the Milky Way (e.g., Mackereth et al., 2018; Gebek & Matthee, 2022). In contrast, the clumpy star formation models predict that chemical bimodality would be more common in galaxies with comparable mass to the Milky Way, because star formation clumps are observed in more than 60% of high-redshift galaxies (e.g., Guo et al., 2015). Such an analysis could be done using spatially-resolved stellar population deconstruction of an edge-on disk, the likes of which are only recently becoming achievable observationally (e.g., Martig et al., 2021).