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

Julie Imig, Kathryn Price, Jon A. Holtzman, Alexander Stone-Martinez, Steven R. Majewski, David H. Weinberg, Jennifer A. Johnson, Carlos Allende Prieto, Rachael L. Beaton, Timothy C. Beers, et al.

► To cite this version:

Julie Imig, Kathryn Price, Jon A. Holtzman, Alexander Stone-Martinez, Steven R. Majewski, et al. A Tale of Two Disks: Mapping the Milky Way with the Final Data Release of APOGEE. The Astrophysical Journal, 2023, 954, 10.3847/1538-4357/ace9b8. insu-04198959

HAL Id: insu-04198959

https://insu.hal.science/insu-04198959

Submitted on 7 Sep 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

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

Julie Imig<fnref n="1" />, Cathryn Price<fnref n="1" />, Jon A. Holtzman<fnref n="1" />, Alexander Stone-Martinez<fnref n="1" />, Steven R. Majewski<fnref n="2" />,
David H. Weinberg<fnref n="3" />,⁴, Jennifer A. Johnson<fnref n="3" />, Carlos Allende Prieto<fnref n="5" />,<fnref n="6" />, Rachael L. Beaton<fnref n="7" />,<fnref n="8" />,<fnref n="9" />, Timothy C. Beers<fnref n="10" />,
Dmitry Bizyaev<fnref n="11" />,<fnref n="12" />, Michael R. Blanton<fnref n="13" />, Joel R. Brownstein¹⁴, Katia Cunha<fnref n="15" />, José G. Fernández-Trincado<fnref n="16" />,
Diane K. Feuillet<fnref n="17" />, Sten Hasselquist<fnref n="7" />, Christian R. Hayes<fnref n="18" />, Henrik Jönsson<fnref n="19" />, Richard R. Lane<fnref n="20" />, Jianhui Lian<fnref n="21" />,
Szabolcs Mészáros<fnref n="22" />,<fnref n="23" />,<fnref n="24" />, David L. Nidever<fnref n="25" />, Annie C. Robin<fnref n="26" />, Matthew Shetrone<fnref n="27" />, Verne Smith<fnref n="28" />, and
John C. Wilson<fnref n="29" />

*Received 2023 May 15; revised 2023 July 21; accepted 2023 July 21; published 2023 August 29*

## Abstract

We present new maps of the Milky Way disk showing the distribution of metallicity ([Fe/H]), α-element abundances ([Mg/Fe]), and stellar age, using a sample of 66,496 red giant stars from the final data release (DR17) of the Apache Point Observatory Galactic Evolution Experiment survey. We measure radial and vertical gradients, quantify the distribution functions for age and metallicity, and explore chemical clock relations across the Milky Way for the low-α disk, high-α disk, and total population independently. The low-α disk exhibits a negative radial metallicity gradient of −0.06 ± 0.001 dex kpc⁻<fnref n="1" />, which flattens with distance from the midplane. The high-α disk shows a flat radial gradient in metallicity and age across nearly all locations of the disk. The age and metallicity distribution functions shift from negatively skewed in the inner Galaxy to positively skewed at large radius. Significant bimodality in the [Mg/Fe]–[Fe/H] plane and in the [Mg/Fe]–age relation persist across the entire disk. The age estimates have typical uncertainties of ~0.15 in log(age) and may be subject to additional systematic errors, which impose limitations on conclusions drawn from this sample. Nevertheless, these results act as critical constraints on galactic evolution models, constraining which physical processes played a dominant role in the formation of the Milky Way disk. We discuss how radial migration predicts many of the observed trends near the solar neighborhood and in the outer disk, but an additional more dramatic evolution history, such as the multi-infall model or a merger event, is needed to explain the chemical and age bimodality elsewhere in the Galaxy.

*Unified Astronomy Thesaurus concepts: Milky Way Galaxy (1054); Milky Way disk (1050); Galactic abundances (2002); Stellar ages (1581); Galaxy stellar content (621); Galactic archaeology (2178); Galaxy structure (622); Milky Way formation (1053); Milky Way evolution (1052)*

<figure><img src="image_3.png" /></figure>

Original content from this work may be used under the terms of the [Creative Commons Attribution 4.0 licence](https://creativecommons.org/licenses/by/4.0/). Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

## 1. Introduction

The positions, chemical compositions, and ages of individual stars in the Milky Way reflect the formation and evolution

<footnote n="1">¹Department of Astronomy, New Mexico State University, P.O. Box 30001, MSC 4500, Las Cruces, NM 88033, USA</footnote>

<footnote n="2">²Department of Astronomy, University of Virginia, Charlottesville, VA 22904, USA</footnote>

<footnote n="3">³The Department of Astronomy and Center of Cosmology and Astro Particle Physics, The Ohio State University, Columbus, OH 43210, USA
⁴Institute for Advanced Study, Princeton, NJ 08540, USA</footnote>

<footnote n="5">⁵Instituto de Astrofísica de Canarias, E-38205 La Laguna, Tenerife, Spain</footnote>

<footnote n="6">⁶Departamento de Astrofísica, Universidad de La Laguna, E-38206 La Laguna, Tenerife, Spain</footnote>

<footnote n="7">⁷Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA</footnote>

<footnote n="8">⁸Department of Astrophysical Sciences, Princeton University, Princeton, NJ 08544, USA</footnote>

<footnote n="9">⁹The Observatories of the Carnegie Institution for Science, 813 Santa Barbara Street, Pasadena, CA 91101, USA</footnote>

<footnote n="10">¹⁰Department of Physics and Astronomy and JINA Center for the Evolution of the Elements, University of Notre Dame, Notre Dame, IN 46556, USA</footnote>

<footnote n="11">¹¹Apache Point Observatory and New Mexico State University, P.O. Box 59, Sunspot, NM 88349-0059, USA</footnote>

<footnote n="12">¹²Sternberg Astronomical Institute, Moscow State University, Moscow</footnote>

<footnote n="13">¹³Center for Cosmology and Particle Physics, Department of Physics, 726 Broadway, Room 1005, New York University, New York, NY 10003, USA
¹⁴Department of Physics and Astronomy, University of Utah, 115 S. 1400 E., Salt Lake City, UT 84112, USA</footnote>

<footnote n="15">¹⁵Steward Observatory, University of Arizona, Tucson, AZ 85721, USA</footnote>

<footnote n="16">¹⁶Instituto de Astronomía, Universidad Católica del Norte, Av. Angamos 0610, Antofagasta, Chile</footnote>

<footnote n="17">¹⁷Lund Observatory, Department of Astronomy and Theoretical Physics, Box 43, SE-221 00 Lund, Sweden</footnote>

<footnote n="18">¹⁸NRC Herzberg Astronomy and Astrophysics Research Centre, 5071 West Saanich Road, Victoria, B.C., V9E 2E7, Canada</footnote>

<footnote n="19">¹⁹Materials Science and Applied Mathematics, Malmö University, SE-205 06 Malmö, Sweden</footnote>

<footnote n="20">²⁰Centro de Investigación en Astronomía, Universidad Bernardo O'Higgins, Avenida Viel 1497, Santiago, Chile</footnote>

<footnote n="21">²¹South-Western Institute for Astronomy Research, Yunnan University, Kunming, Yunnan 650091, People's Republic of China</footnote>

<footnote n="22">²²ELTE Eötvös Loránd University, Gothard Astrophysical Observatory, 9700 Szombathely, Szent Imre H. st. 112, Hungary</footnote>

<footnote n="23">²³MTA-ELTE Lendület "Momentum" Milky Way Research Group, Hungary</footnote>

<footnote n="24">²⁴MTA-ELTE Exoplanet Research Group, Szombathely, Szent Imre h.u. 112., H-9700, Hungary</footnote>

<footnote n="25">²⁵Department of Physics, Montana State University, P.O. Box 173840, Bozeman, MT 59717, USA</footnote>

<footnote n="26">²⁶Institut UTINAM—UMR 6213—CNRS—University of Bourgogne Franche Comté, France, OSU THETA, 41bis avenue de l'Observatoire, F-25000, Besançon, France</footnote>

<footnote n="27">²⁷University of California Observatories, UC Santa Cruz, Santa Cruz, CA 95064, USA</footnote>

<footnote n="28">²⁸NSF's National Optical-Infrared Astronomy Research Laboratory, 950 North Cherry Avenue, Tucson, AZ 85719, USA</footnote>

<footnote n="29">²⁹Astronomy Department, University of Virginia, Charlottesville, VA 22901, USA</footnote>
history of our Galaxy, with each individual star acting as a “fossil” containing the chemical fingerprint of the interstellar gas from which it formed. Our inside perspective in the Milky Way grants us the ability to study it in greater detail than any other galaxy, placing strong observational constraints on formation models and simulations of disk galaxies. For this reason, constraining the chemical and dynamical properties of the stellar populations in the Milky Way disk remains a cornerstone of modern galactic astronomy.

Understanding the present-day chemical structure of our Galaxy has been increasingly successful with the advent of large spectroscopic stellar surveys like Apache Point Observatory Galactic Evolution Experiment (APOGEE; Majewski et al. 2017), Gaia (Gaia Collaboration et al. 2021), Gaia-ESO (Gilmore et al. 2022), LAMOST (Luo et al. 2015), GALAH (Buder et al. 2018), RAVE (Steinmetz et al. 2020), and SEGUE (e.g., Yanny et al. 2009). These surveys obtain precise kinematic and chemical information for a combined millions of stars across the Milky Way, with increasing sample sizes and more complete spatial coverage with every generation of survey. When paired with precise distances and positions from Gaia astrometry (Gaia Collaboration et al. 2018, 2021), these large surveys can access the evolution history of a large fraction of the Galactic disk. Even stellar ages, notoriously difficult to infer as they previously could not be directly measured for individual stars, are now readily available through the precise measurements of the masses for thousands of red giant stars available through asteroseismology (e.g., Pinsonneault et al. 2018; Miglio et al. 2021). These asteroseismic data sets are additionally used as training sets for machine-learning techniques, expanding the stellar sample with age estimates to hundreds of thousands of stars (e.g., Ness et al. 2016; Anders et al. 2018; Leung & Bovy 2018; Mackereth et al. 2019; Wu et al. 2019; Ciucă et al. 2021; Stone-Martinez et al. 2023, submitted).

Despite this wealth of data, the debate remains heated around which physical processes played the largest roles in shaping the Milky Way's disk. The structural and chemical distribution of stars in the Milky Way has been well studied, leading to the discovery of two main stellar components, the “thin” and “thick” disk near the solar neighborhood (e.g., Yoshii 1982; Gilmore & Reid 1983). These components are distinct in their dynamic signature, with the thick disk characterized by kinematically hotter stellar orbits (larger vertical velocity dispersion), and a slower systemic rotational velocity than the thin disk (e.g., Soubiran et al. 2003; Jurić & Ivezić et al. 2008; Kordopatis et al. 2013; Robin et al. 2017). The thin disk is also generally accepted to be more radially extended, and as the name implies, has a smaller scale height than the thick disk (e.g., Bensby et al. 2011; Bovy et al. 2016; Mackereth et al. 2017; Lian et al. 2022b; Robin et al. 2022). The two disks also differ in their chemical fingerprints, with the thin disk generally containing younger metal-rich stars characterized by their lower α-element<fnref n="30" /> abundances relative to the older, more metal-poor thick disk (e.g., Fuhrmann 1998; Bensby et al. 2005; Reddy et al. 2006; Lee et al. 2011a; Bovy et al. 2012b, 2016; Nidever et al. 2014; Hayden et al. 2015; Kordopatis et al. 2015a; Mackereth et al. 2017; Katz 2021; Vincenzo et al. 2021a). The [Mg/Fe] ratio reflects the relative iron enrichment by prompt, massive core-collapse supernovae compared to the longer timescale Type Ia supernovae. Because

of this, the [Mg/Fe] ratio is generally high in populations that formed during rapid and efficient starbursts, and approaches solar “α-poor” values in populations that form steadily over long time periods (e.g., Matteucci & Brocato 1990; Thomas et al. 2005). Thus, the chemical differences between the thin and thick disk suggests that they formed via distinct pathways, leaving the evidence of their enrichment histories within the present-day chemical structure of the Galaxy. However, many of these studies are biased toward the solar neighborhood due to observation limitations, and there has been some debate on whether the two components are truly distinct at all (e.g., Bensby et al. 2007; Bovy et al. 2012a; Kawata & Chiappini 2016; Hayden et al. 2017; Anders et al. 2018).

Nevertheless, different explanations for the origins of this chemical bimodality in the disk have been proposed, using a combination of physical processes such as star formation, gas accretion, quenching, galaxy mergers, and stellar radial migration to attempt to explain the observed trends. The different models can be generally categorized into three scenarios: the “two-infall” models where the thick disk forms first followed by the thin disk, the “superposition” models where the two disks form in different parts of the Galaxy and mix through stellar migration, and the “clumpy formation” models where the two disks form simultaneously but with different star formation efficiencies.

The “two-infall” class of models, originally of Chiappini et al. (1997) and Chiappini et al. (2001), describe a scenario wherein the Milky Way first forms from the collapse of primordial gas, creating the progenitor of the present-day thick disk in a fast burst of star formation. The gas reservoir of the Galaxy is then quenched, entering a quiescent period of little star formation until the Galaxy receives a second infall of pristine gas. This accretion of fresh material dilutes the metallicity of the interstellar medium (ISM) before reigniting star formation that forms the thin disk. The second gas infall happens over a longer timescale, allowing for a period of more continuous star formation, resulting in the α-poor nature of thin disk. Linden et al. (2017) constrained the timing of the second infall to be between 7 and 8 Gyr ago based on the ages and chemistry of star clusters in APOGEE. Spitoni et al. (2019a, 2020, 2021) expanded on this model, constraining the length of the delay between the two episodes of gas infall to be between 3 and 5.5 Gyr, and proposing the second gas infall corresponds to a merger event with a gas-rich dwarf galaxy around 8–11 Gyr ago. This may coincide with the Milky Way’s accretion of the Gaia-Enceladus dwarf galaxy, estimated to have happened 10 Gyr ago (Helmi et al. 2018; Vincenzo et al. 2019).

A number of three-infall models have also recently been proposed, including the model of Spitoni et al. (2022b) constrained to Gaia data. Their most recent infall starts ~2.7 Gyr ago and gives birth to the recently discovered young, low-α stars that are impoverished in some elements (Gaia Collaboration et al. 2023). This latest infall may be linked with the Sagittarius dwarf spheroidal galaxy’s most recent perigalactic passage through the Milky Way’s disk (Laporte et al. 2019; Antoja et al. 2020; Ruiz-Lara et al. 2020). A starburst 2–3 Gyr ago has been detected independently in Isern (2019) and Mor et al. (2019).

The works of Lian et al. (2020a, 2020b, 2020c) and Lian et al. (2021) present a modified version of the two-infall model. In their version, an underlying continuous episode of gas accretion is interrupted by two rapidly quenched starbursts. The first starburst forms the high-α thick disk, and the second

<footnote n="30">³⁰ α-elements are elements with an atomic number multiple of 4 (the mass of a Helium nucleus, an α-particle), e.g., O, Mg, S, and Ca.</footnote>
starburst forms the metal-poor end of the low-$\alpha$ sequence 6 Gyr later. The metal-rich low-$\alpha$ sequence is attributed to the secular evolution phase between the two bursts.

Another variation of the two-infall model without the inclusion of merger events has been supported by recent chemo-dynamical simulations from Khoperskov et al. (2021). As in previous models, the thick disk is formed early on in a burst of star formation in a turbulent, compact disk. Stellar feedback from the formation of the thick disk drives outflows that quench star formation, enrich the Galactic halo, and eventually feed the gas back into the disk on a more sustained timescale, creating the thin disk with a “galactic fountain” (e.g., Shapiro & Field 1976; Bregman 1980; Marinacci et al. 2011; Fraternali 2017). The models of Haywood et al. (2016), Haywood et al. (2018), and Haywood et al. (2019) support this scenario, where the high-$\alpha$ population was formed early on in a turbulent gas-rich disk with strong feedback, and the leftover, diluted gas forms the low-$\alpha$ thin disk on longer timescales.

The “superposition” class of chemical evolution models, pioneered by Schönrich & Binney (2009a) and Schönrich & Binney (2009b), reproduces the observed disk dichotomy without the need for a violent merger history to heat the thick disk. In this scenario, the chemical locus of the thin disk is not an evolutionary track; it is a superposition of end points of evolutionary tracks from different Galactocentric radii (e.g., Nidever et al. 2014; Kubryk et al. 2015; Sharma et al. 2021a). Stars from these different tracks reach the solar neighborhood by radial migration, a natural consequence of the Galaxy’s spiral structure (Sellwood & Binney 2002; Roškar et al. 2008). Stars in the high-$\alpha$ thick disk formed early during an efficient phase of rapid star formation, primarily in the inner Galaxy, before migrating to their present-day radial distribution.

This superposition model is expanded on in the works of Minchev et al. (2013, 2014), Minchev et al. (2017), and Johnson et al. (2021), which also emphasize the importance of radial migration in the Milky Way’s structure. In these models, gas inflows, outflows, and star formation rates vary with Galactic location, emphasizing the radial dependence of the disk’s chemical evolution history. Stellar radial migration allows stars to move around the Galaxy as time progresses, and potentially enrich a different spatial zone than the one they were born in when they die. These models show that this radial migration is the key to reproducing many of the observed trends in the Milky Way, including the changes of [Mg/Fe] and [Fe/H] distributions with radius and height.

A third, qualitatively different scenario is proposed by the “clumpy formation” model of Clarke et al. (2019), motivated by results from hydrodynamical simulations (e.g., Bournaud et al. 2007) and observations of high-redshift galaxies (e.g., Elmegreen et al. 2005). In their picture, the low-$\alpha$ thin disk is a true evolutionary sequence corresponding to inefficient star formation in the disk, while the high-$\alpha$ population is formed mainly during rapid, clumpy bursts in the Galaxy’s early gas-rich phase. These clumps are comparable to those observed in high-redshift galaxies with the Hubble Space Telescope (e.g., Elmegreen et al. 2005). In addition to the chemical bimodality, these models also reproduce the observed mass density structure of the Milky way, including the flared thin disk (Amarante 2020; Beraldo e Silva et al. 2020).

The discrepancy between these different proposed explanations, which all reasonably reproduce the observed trends in the Milky Way’s disk, can only be closed with more observational

constraints. Detailed chemical maps that cover the entire span of the disk, robust measurements of the Milky Way’s radial and vertical metallicity gradients, and the metallicity distribution function (MDF) will help constrain which physical processes played an important role in the formation of the disk. Adding in the ages of stars can provide an important axis for interpreting these results, as they enable a direct temporal connection between the properties of individual stars and the evolutionary timescale of the Milky Way (e.g., Mackereth et al. 2017; Feuillet et al. 2019; Vázquez et al. 2022).

In this paper, the final data release of APOGEE (Majewski et al. 2017; Abdurro’uf et al. 2022) is used to further explore the properties of the Milky Way, with a larger sample size and greater spatial coverage than previously available.

The chemical cartography of the Milky Way has been extensively studied previously using a variety of different surveys including SEGUE (e.g., Lee et al. 2011a, 2011b; Gómez et al. 2012), RAVE (e.g., Kordopatis et al. 2013; Robin et al. 2017), GALAH (e.g., Lin et al. 2019; Hayden et al. 2020; Sharma et al. 2021b), LAMOST (e.g., Huang et al. 2020; Vickers et al. 2021; Hawkins 2023), Gaia (e.g., Lemasle et al. 2018; Poggio et al. 2022; Gaia Collaboration et al. 2023), Gaia-ESO (e.g., Bergemann et al. 2014; Kordopatis et al. 2015a; Magrini et al. 2018; Vázquez et al. 2022), and previous data releases of APOGEE (e.g., Nidever et al. 2014; Hayden et al. 2015; Weinberg et al. 2019; Katz 2021; Eilers et al. 2022). These works, and others, have impressively advanced the field of chemical cartography over the last decade, meaning that many of the results presented in this paper are not new. However, the final data release of APOGEE presents a larger and more detailed base data set than previously available. This allows us to cull a selected high-quality sample, minimizing systematics while still probing a large number of stars at different locations across the Galaxy. Additionally, APOGEE has the distinct advantage of working in the infrared, easily accessing the heavily dust-obscured regions like the Galactic center and midplane, which are often beyond the reach of optical surveys.

Our study complements the DR17-based study of Weinberg et al. (2022), which focused on abundance trends of [X/Mg] for many different elements. These trends, which are nearly universal throughout the disk, provide insights on nucleosynthetic processes, while the distribution of stars in [Mg/Fe], [Fe/H], and age across the disk provide constraints on Galactic history.

In this work, we build upon the decades of previous discoveries and explore the chemical trends in the Milky Way disk through the legacy of the APOGEE survey. A high-quality sample of 66,496 red giant stars and their precise measurements of metallicity ([Fe/H]), ages, and $\alpha$-element abundances ([Mg/Fe]) are used to create maps, measure gradients, quantify distribution functions, and trace age-abundance relations across the Milky Way disk, and compare the observations with the most recent models. In short, we find evidence supporting all three classes of chemical evolution models. Radial migration is an important process in shaping the disk over time, but the observed bimodality in $\alpha$-element abundances and ages persists even in disk regions where radial migration is not expected to be as prevalent. This suggests a multiphase star formation history, such as that presented in the two-infall or clumpy formation class of models, is at least
partially responsible for the formation of the Milky Way as seen today.

Section 2 contains an overview of the APOGEE survey and supplementary data used in this study. Spatial maps, gradient measurements, distribution functions, and other results are presented in Section 3 and compared with previous literature. In Section 4, we discuss our results in the context of chemical evolution models. The conclusions we draw from this study are presented in Section 5.

## 2. Data

### 2.1. APOGEE

APOGEE (Majewski et al. 2017) is a high-resolution ($R \sim 22,500$) near-infrared (1.51–1.70 µm) spectroscopic survey containing observations of 657,135 unique stars released as part of the Sloan Digital Sky Survey (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 2021 December as part of the 17th data release (DR17) of SDSS (Abdurro'uf et al. 2022) and is available publicly online through the SDSS Science Archive Server and Catalog Archive Server.<fnref n="31" />

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 preparation). ASPCAP derives stellar atmospheric parameters, radial velocities, and as many as 20 individual elemental abundances for each APOGEE spectrum by comparing each to a multidimensional 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 α-element abundances ([Mg/Fe]) from ASPCAP. We adopt [Mg/Fe] for our α-element abundance instead of the “total” [α/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.

### 2.2. Sample Selection

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 Two Micron All Sky Survey 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 \le \log g \le 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_{\text{eff}}$ and $\log g$ in APOGEE (e.g., Jönsson et al. 2018; Eilers et al. 2022). 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. (2022) presented 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_{\text{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 \ge 0; N = 38,031$) than below ($Z \le 0; N = 28,465$), and more observations toward the Galactic center ($R \le R_{\odot}; N = 36,317$) than outward ($R \ge 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 remainder 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.<fnref n="32" /> 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 toward the solar neighborhood. Additionally, nothing can be inferred from the

<footnote n="31">³¹ Data Access Instructions: https://www.sdss.org/dr17/irspec/spectro_data/.</footnote>

<footnote n="32">³² This 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 nonhomogeneous dust distribution in the Milky Way plays a major role in what can be observed.</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 1.** The $T_{eff}$-log g distribution of stars in the sample described in Section 2. Our adopted red giant sample is outlined in black and plotted by color (metallicity), while the full APOGEE sample is shown in gray in the background for reference.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 2.** The spatial distribution of stars in our sample shown as a face-on view of the Galaxy (X-Y plane; top panel), and an edge-on view (R-Z plane; bottom panel). Each spatial bin of $\Delta X = \Delta Y = \Delta Z=0.5$ kpc is colored by the number of stars in that location. The position of the Sun at $X = -8.3$ kpc is denoted by the solar symbol (☉), and the Galactic center is marked with a plus (+). The ellipse around the Galactic center marks the approximate location of the bar, as an ellipse with major axis length 10 kpc, a 0.4 axis ratio, and rotated 25°.</figcaption></figure>

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

<figure><img src="image_6.png" /><figcaption>**Figure 3.** The [Mg/Fe]-[Fe/H] plane for stars in our red giant sample (black points), demonstrating the adopted definition of the low-$\alpha$ (blue) and high-$\alpha$ (red) sequence defined in Equation (1). The gray region is an added buffer zone of [Fe/H] = ±0.025 dex around the line to remove overlap between the two sequences due to abundance uncertainties. The typical uncertainties associated with each measurement of [Fe/H] and [Mg/Fe] are shown in the error bars near the bottom of the plot (as the ±1σ value), for several selected metallicities.</figcaption></figure>

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) and Mackereth et al. (2017) for previous data releases, and J. Imig et al. (2023, in preparation) for DR17. J. Imig et al. (2023, in preparation) will present the density distribution of mono-age mono-abundance stellar populations in APOGEE DR17 after correcting for the selection function.

## 2.3. [Mg/Fe] Subsamples

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. 2011a; Kordopatis et al. 2015a; Hayden et al. 2015; Katz 2021). Although there is some debate on whether the two sequences are truly distinct (see Section 1), we use this figure to define two further subsamples in our data to investigate this question later. Vincenzo et al. (2021a) demonstrated 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:

$$ \text{Mg/Fe} = \begin{cases} 0.15 - 0.13*[\text{Fe/H}] & \text{if } [\text{Fe/H}] \le 0 \\ 0.15 & \text{if } [\text{Fe/H}] > 0 \end{cases} \quad (1) $$

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 ±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 ±1σ in the error bars across the bottom of the plot.

## 2.4. Age and Distance Estimates

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 R.A. and decl. from APOGEE observations and distance estimates from the APOGEE *distmass* value added catalog (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_⊙ = 8.3$ kpc with a height of $z_⊙ = 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: M.H. Pinsonneault et al. 2023, in preparation). 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 \le \log(\text{age}) \le 10.3$), metallicities ($-2.0 \le [\text{Z/H}] \le 0.5$), and masses ($0.1 \le M/M_⊙ \le 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}(\text{age})$, and the median upper uncertainty is $0.16 \log_{10}(\text{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$ (Stone-Martinez et al. 2023, submitted).

To determine if these uncertainties are realistic, Stone-Martinez et al. (2023, submitted) performed 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 found that the *distmass* ages are accurate within ±1σ = 0.16 in $\log_{10}(\text{age})$ across 12 different star clusters with ages $9.2 \le \log_{10}(\text{age}) \le 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 ±1σ = 0.11 in $\log_{10}(\text{age})$, although these ages are derived with a similar methodology. Both of these evaluations are consistent within the reported age uncertainties.

<figure><img src="image_7.png" /><figcaption>**Figure 4.** The reported ±1σ uncertainty in stellar age estimates from the *distmass* value added catalog. The light-gray histogram shows the lower uncertainty values for the sample, and the dark-gray histogram shows the upper uncertainties for the sample. The median of each distribution is plotted as a vertical line and labeled to highlight a “typical” uncertainty value.</figcaption></figure>

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. The Appendix 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 (M.H. Pinsonneault et al. 2023, in preparation), 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 (M.H. Pinsonneault et al. 2023, in preparation), 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 the Appendix. 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 the Appendix.

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_{\text{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 $[\text{Fe/H}] \le -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.
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Global maps of the Milky Way, showing the average distribution of [Fe/H] (left), [Mg/Fe] (middle), and stellar age (right) across the Galaxy. The top row of panels shows a face-on view (X-Y plane), integrated through the disk with $|Z| \le 1.0$ kpc. The bottom row of panels shows an edge-on view (R-Z plane, with R preserving the sign of X to show the opposite side of the Galaxy), integrated through the whole disk. Colors encode the median quantities in each or X-Y or R-Z pixel. In the face-on views, the age and metallicity gradients are visible, with the Galactic bar standing out as metal-rich and α-poor. The location of the Sun at X = −8.3 kpc is denoted by the solar symbol (☉), and the Galactic center is marked with a plus (+). The approximate location of the Galactic bar is also shown as an ellipse with major axis length 10 kpc, a 0.4 axis ratio, and rotated 25°.</figcaption></figure>

Because metal-poor stars tend to be older, this means that for age-related figures, the oldest stars (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 the Appendix.

## 3. Results

### 3.1. Cartography

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 toward outer radii. In median α-element abundances (middle column), the bimodality in the disk shows low-α stars congregating in the “thin disk” near the Galactic midplane, and high-α stars populating the “thick disk” at higher Z locations. At larger radii, the low-α stars extend farther above and below the plane. The high-α 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, α-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; Lian et al. 2021; Queiroz et al. 2021; Eilers et al. 2022) 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| \le 0.5$ kpc) to farthest away (bottom panel; $1 \le |Z| \le 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 midplane, the trend becomes less apparent, with almost no obvious gradient present when $|Z| \ge 1$ kpc, and the stellar populations showing a lower median metallicity overall. The middle panel ($0.5 \le |Z| \le 1$ kpc) shows a peculiar trend where the median metallicity actually increases from R=0 until R ~ 7 kpc, and then decreases with a shallow metallicity gradient. This buildup of metal-rich stars in the center of the Galaxy is possibly a 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.

### 3.2. [Mg/Fe] Distribution

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
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Face-on maps of the Galactic disk, showing the median metallicity ([Fe/H]) of stars in spatial bins of $\Delta X = \Delta Y = 0.5$ kpc. The different panels are slices in vertical space, from closest to the Galactic plane (top panel; $|Z| < 0.5$ kpc) to increasing heights above the plane (bottom panel; $1 < |Z| < 2$ kpc). The Sun's position is marked by the solar symbol (☉), and the position of the Galactic center is indicated by a plus (+).</figcaption></figure>

(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.

<figure><img src="image_10.png" /><figcaption>**Figure 7.** Same as Figure 6, but colored by the median age (in gigayears) of stars in each spatial bins of $\Delta X = \Delta Y = 0.5$ kpc.</figcaption></figure>

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
<figure><img src="image_11.png" /><figcaption>**Figure 8.** The distribution of stars in the [Mg/Fe] vs. [Fe/H] plane as a function of R and |Z|, as a contour map of point density. Spatial bins move from closest to the Galactic plane (bottom row, 0.0 < |Z| < 0.5 kpc) to farthest above the Galactic plane (top row, 1.0 < |Z| < 2.0 kpc), and from close to the Galactic center (left column, 0.0 < |R| < 3.0 kpc) to farthest out in the disk (right column, 15.0 < |R| < 25.0 kpc). The number in the top-right corner of each panel is the number of stars in our sample in that spatial bin. For reference, the gray background shape and black line are the same in each panel, to highlight how the sequence changes across location in the Galaxy. The black line is the boundary between high- and low-α populations defined in Equation (1), and the gray shape is the contour containing 90% of the points in the full sample. The typical uncertainties in abundance measurements as a function of metallicity are shown as a ±1σ value at the bottom of each panel for reference.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 9.** Same as Figure 8, but colored by stellar age. Points with [Fe/H] ≤ −0.7 have been excluded for potentially unreliable age estimates with the cuts described in Section 2.4. The percentage of stars ($f_{distmass}$) in each bin that pass the distmass quality criterion is shown in parentheses in the upper-right corner of each panel.</figcaption></figure>

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 2021; Vincenzo et al. 2021a; Gaia Collaboration et al. 2023).

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 outward 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 farther 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 downward (toward more $\alpha$-poor) with increasing radius. This is also seen in Katz (2021) with APOGEE data, using the mode of the data, although they found a larger shift of ~0.05 dex between the inner and outer Galaxy, while ours is closer to half that at ~0.025 dex.

### 3.3. Azimuthal Variance in Metallicity

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 nonaxisymmetric features such as the Galactic
<figure><img src="image_13.png" /><figcaption>**Figure 10.** Variation in the [Mg/Fe]–[Fe/H] plane as a function of Galactic radius (color) for stars within $|Z| < 0.5$ kpc of the plane. Top panel: the contour containing 75% of points for both the low-$\alpha$ and high-$\alpha$ samples, calculated separately. Bottom panel: the median [Mg/Fe] as a function of [Fe/H] for both sequences. The shaded regions denotes the $\pm 1\sigma$ uncertainty in the median. In both panels, the black line is our defined boundary between the low-$\alpha$ and high-$\alpha$ described in Equation (1).</figcaption></figure>

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^\circ \le \theta \le 230^\circ$, 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_{[Fe/H]}$) is defined for a given radius as the sum of the uncertainty in the median for each individual $\theta$ bin ($\sigma_{[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:

$$ \sigma_{[Fe/H]}(R)^2 = \frac{1}{N_{\text{bins},\theta}} \sum_{\theta=0}^{360} \sigma_{[Fe/H]}(R, \theta)^2 \quad (2) $$

$$ \sigma_{[Fe/H]}(R, \theta)^2 = \frac{1}{N_{\text{stars}}} \sum_{i=0}^{N_{\text{stars}}} |([\text{Fe/H}]_i - \overline{[\text{Fe/H}]})|^2. \quad (3) $$

Close to the Galactic plane ($|Z| \le 0.5$ kpc, top row), the Galaxy has little metallicity variation across azimuth near the

solar neighborhood and outward ($R \ge 8$ kpc), although the region of azimuth covered by the observations decreases with increasing radius. Closer to the center of the Galaxy ($R \le 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 \le |Z| \le 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 \le |Z| \le 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 nonaxisymmetric ways, such as warping the disk or introducing kinematic oscillations (e.g., Gómez et al. 2012; Cheng et al. 2020; Chrobáková 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 2023) possibly corresponding to the Milky Way’s spiral arms. Measuring the radial metallicity gradient using Gaia data, Poggio et al. (2022) and Hawkins (2023) reported 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 a 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.

### 3.4. Radial and Vertical Metallicity Gradients

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; Carrell et al. 2012; Cheng et al. 2012; Anders et al. 2014; Hayden et al. 2014; Schlesinger et al. 2014; Frankel et al. 2019; Katz 2021; Gaia Collaboration et al. 2023). Such a trend is predicted by “inside-out” disk formation models, where stars
<figure><img src="image_14.png" /><figcaption>**Figure 11.** Maps of the median metallicity in the disk in polar coordinates, to highlight any nonaxisymmetric features in the disk. The spatial bins are sized $\Delta R = 1$ kpc and $\Delta\theta = 10^{\circ}$. The rows are different slices in $Z$, moving from closest to the Galactic plane (top row, $|Z| \le 0.5$ kpc), to farther above (bottom row, $1 \le |Z| \le 2.0$ kpc). The left column shows the median [Fe/H] for each bin as a function of $\theta$ (point color), compared to the $\pm 1\sigma$ expected spread from uncertainty in the measurement of the median (gray shaded region). The right column shows the median metallicity maps as a function of radius $R$ and azimuthal angle $\theta$. The black line is the approximate location of the Galactic Bar, defined as an ellipse with major axis length 10 kpc, a 0.4 axis ratio, and rotated $25^{\circ}$. The black square in the top panel highlights the region shown in more detail in Figure 12.</figcaption></figure>

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. 2022a). 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á & Díaz 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).
<figure><img src="image_15.png" /><figcaption>**Figure 12.** Azimuthal variation in metallicity restricted to a window near the solar neighborhood, outlined by the black rectangle in Figure 11. The metallicity value is calculated as a running median for bin size $\Delta R = 0.5$ kpc and $\Delta\theta = 5^{\circ}$, evaluated every $\Delta R = 0.1$ kpc and $\Delta\theta = 1^{\circ}$ for smoothing. The approximate locations of the nearby spiral arms from Reid et al. (2019) are plotted as colored lines.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Figure 13.** Radial median metallicity profile as a function of height out of the plane (line color), for the total stellar population (top panel), the low-$\alpha$ disk (middle panel), and the high-$\alpha$ disk (bottom panel).</figcaption></figure>

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.

<table><caption>**Table 1**

Radial Metallicity Gradients in dex kpc$^{-1}$ as a Function of Height above the Plane |Z| from Figure 14</caption><thead><tr><th>|Z| (kpc)</th><th>Total</th><th>Low-α</th><th>High-α</th></tr></thead><tbody><tr><td>0.0</td><td>-0.056 ± 0.001</td><td>-0.06 ± 0.001</td><td>-0.01 ± 0.007</td></tr><tr><td>0.25</td><td>-0.057 ± 0.001</td><td>-0.064 ± 0.001</td><td>-0.008 ± 0.006</td></tr><tr><td>0.5</td><td>-0.051 ± 0.001</td><td>-0.063 ± 0.001</td><td>-0.01 ± 0.006</td></tr><tr><td>0.75</td><td>-0.041 ± 0.002</td><td>-0.054 ± 0.002</td><td>-0.006 ± 0.005</td></tr><tr><td>1.0</td><td>-0.023 ± 0.002</td><td>-0.04 ± 0.002</td><td>-0.0 ± 0.005</td></tr><tr><td>1.25</td><td>-0.01 ± 0.002</td><td>-0.023 ± 0.003</td><td>0.003 ± 0.006</td></tr><tr><td>1.5</td><td>-0.002 ± 0.003</td><td>-0.005 ± 0.005</td><td>0.006 ± 0.006</td></tr><tr><td>1.75</td><td>0.001 ± 0.003</td><td>0.013 ± 0.005</td><td>-0.01 ± 0.007</td></tr></tbody></table>

The total sample (Figure 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 (Figure 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 (Figure 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 \ge 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 \ge 7$ kpc), the high-$\alpha$ population shows 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| ≤ 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. (2023) 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)
<figure><img src="image_17.png" /><figcaption>**Figure 14.** The best-fit slope for each radial metallicity profile in Figure 13 in units of dex kpc⁻¹, fit with a single line for stars beyond R > 7 kpc for the total sample (black line), and the low-α (blue line) and high-α (red line) samples independently. The shaded region indicates the ±1σ uncertainty in the slope measurement.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Figure 15.** Vertical median metallicity profile as a function of Galactocentric radius (line color) for the total stellar population (top panel), the low-α disk (middle panel), and the high-α disk (bottom panel).</figcaption></figure>

reporting −0.060 ± 0.002 dex kpc⁻¹, and Lemasse et al. (2018) reporting −0.045 ± 0.007 dex kpc⁻¹. The large sample size, distance range, and high precision of the APOGEE DR17 sample enable us to map the radial, vertical, and α-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 midplane, which alone results in a steeper vertical gradient.

<table><caption>**Table 2**
Vertical Metallicity Gradients in dex kpc⁻¹ as a Function of Galactocentric Radius *R* from Figure 16</caption><thead><tr><th><i>R</i> (kpc)</th><th>Total</th><th>Low-α</th><th>High-α</th></tr></thead><tbody><tr><td>0.0</td><td>−0.4 ± 0.053</td><td>−0.565 ± 0.063</td><td>−0.126 ± 0.042</td></tr><tr><td>2.0</td><td>−0.471 ± 0.015</td><td>−0.616 ± 0.027</td><td>−0.121 ± 0.012</td></tr><tr><td>4.0</td><td>−0.462 ± 0.01</td><td>−0.753 ± 0.021</td><td>−0.084 ± 0.009</td></tr><tr><td>6.0</td><td>−0.444 ± 0.011</td><td>−0.519 ± 0.019</td><td>−0.1 ± 0.013</td></tr><tr><td>8.0</td><td>−0.296 ± 0.009</td><td>−0.265 ± 0.013</td><td>−0.119 ± 0.013</td></tr><tr><td>10.0</td><td>−0.153 ± 0.007</td><td>−0.145 ± 0.008</td><td>−0.062 ± 0.015</td></tr><tr><td>12.0</td><td>−0.09 ± 0.007</td><td>−0.094 ± 0.007</td><td>−0.075 ± 0.027</td></tr><tr><td>14.0</td><td>−0.05 ± 0.009</td><td>−0.049 ± 0.009</td><td>−0.115 ± 0.071</td></tr><tr><td>16.0</td><td>−0.066 ± 0.017</td><td>−0.067 ± 0.017</td><td>−0.033 ± 0.792</td></tr><tr><td>18.0</td><td>−0.072 ± 0.04</td><td>−0.072 ± 0.041</td><td></td></tr></tbody></table>

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-α population (Hayden et al. 2014, their 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 ± 0.001 dex kpc⁻¹ is slightly shallower than the −0.073 ± 0.003 dex kpc⁻¹ 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. 2022b; Robin et al. 2022; Gaia Collaboration et al. 2023).

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. (2023), 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-α 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-α population has a shallow negative gradient everywhere, which does not significantly change with radius. Beyond R ≳ 10 kpc, the vertical gradients for the total, high-α, and low-α populations are all close to zero.

This is also consistent with previous studies (e.g., Hayden et al. 2014, their Table 1), where the vertical metallicity
<figure><img src="image_19.png" /><figcaption>**Figure 16.** The best-fit slope for each vertical metallicity profile in Figure 15, fit with a single line for stars beyond $R > 7$ kpc for the total sample (black line), and the low-$\alpha$ (blue line) and high-$\alpha$ (red line) samples independently. The shaded region indicates the $\pm 1\sigma$ uncertainty in the slope measurement.</figcaption></figure>

gradient approaches 0 as radius increases. Near the solar neighborhood, we report a vertical median metallicity gradient of $-0.315 \pm 0.009 \text{ dex kpc}^{-1}$. Hayden et al. (2014) reported a gradient of $-0.31 \pm 0.01 \text{ 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 \le |Z| \le 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 ISM is expected to flatten out over time (e.g., Pilkington & Gibson 2012; Gibson et al. 2013; Mollá & Díaz 2018). The opposite trend, seen here, is commonly attributed as 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.

### 3.5. Metallicity Distribution Function

While the radial and vertical median metallicity gradients in the disk reveal interesting general trends, more insights can be gleaned from the full 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 MDF varies with radius for samples at different vertical slices, from closest to the Galactic plane (top panel) to the farthest 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

<figure><img src="image_20.png" /><figcaption>**Figure 17.** Top panel: the radial $R$ median metallicity profile of the low-$\alpha$ disk, split into different samples of stellar age (line color). Bottom panel: the vertical ($Z$) median metallicity profile of the low-$\alpha$ disk, split into different samples of stellar age (line color).</figcaption></figure>

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 toward lower metallicities (the peak trends right of the median), and in the outer Galaxy, the distribution is skewed toward 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 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,
<figure><img src="image_21.png" /><figcaption>**Figure 18.** The metallicity distribution function (MDF) of the Milky Way disk, split into different height bins (panels), from closest to the Galactic plane (top) to the farthest beyond (bottom). Each panel shows the fraction of stars at each metallicity [Fe/H] as a function of Galactocentric radius, further split by color between high-$\alpha$ (red) and low-$\alpha$ (blue) samples. Every third row is annotated with markings for the peak (or mode) of the distribution (white diamond), as well as the 25th, 50th (median), and 75th percentiles (white tick marks) to highlight the shape of the distribution.</figcaption></figure>

this could imply that the high-$\alpha$ sample is more “mixed” vertically, meaning the birth locations of stars tend to be farther away from their present-day locations, 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 linked it to radial migration. Hayden et al. (2015) observed this same trend in APOGEE and presented 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) showed 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 is negatively skewed, and in the outer regions, it is positively skewed, quantifying the trend seen earlier in Figure 18.

### 3.6. Age Gradients

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 \ge 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, ~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] $\ge -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] $\ge -0.7$ and are therefore more robust against potential biases caused by the sample selection (see Figure A1).

### 3.7. Age Distribution Function

As before with metallicities, more information lies in the shape of the age distribution function (ADF) at different locations in the Galaxy rather than the gradient alone. Figure 22 depicts the ADF
<figure><img src="image_22.png" /><figcaption>**Figure 19.** The MDF and its first three moments, limited to the low-$\alpha$ sample and close to the Galactic plane ($|Z| < 0.5$ kpc; equivalent to the top panel in Figure 18). Left panel: the MDF at different radii in the Galaxy (colored lines). The right panels show the first three statistical moments for quantifying this distribution as a function of radius. Top right: a characteristic [Fe/H] for each distribution function, measured as a median (blue line), a mean (green line), and the peak (or mode; purple line). Middle right: the width, $\sigma$, of the distribution. Bottom right: the skewness of the distribution, where a negative number indicates a left-leaning distribution as shown in the left panel, and a positive number corresponds to a right-leaning distribution. The solar position ($R = 8.3$ kpc) is marked by a vertical gray line in all right panels.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>**Figure 20.** Radial median age profile as a function of height above the plane (line color), for the total stellar population (top panel), the low-$\alpha$ disk (middle panel), and the high-$\alpha$ disk (bottom panel).</figcaption></figure>

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

<caption>**Figure 21.** The best-fit slope for each radial age profile in Figure 20 in units of Gyr kpc$^{-1}$, fit with a single line for stars beyond $R > 7$ kpc for the total sample (black line), and the low-$\alpha$ (blue line) and high-$\alpha$ (red line) samples independently. The shaded region indicates the $\pm 1\sigma$ uncertainty in the slope measurement.</caption>

radii. In the outer Galaxy, the ADF is preferentially skewed toward older ages. Katz (2021) found the ADF skewed toward younger ages in the inner Galaxy, and skewed toward 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 Figure 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
<figure><img src="image_24.png" /><figcaption>**Figure 22.** The age distribution function (ADF) of the Milky Way disk, split into different height bins (panels), from closest to the Galactic plane (top) to farthest beyond (bottom). Each panel shows the fraction of stars at each age as a function of Galactocentric radius, further split by color between high-$\alpha$ (red) and low-$\alpha$ (blue) samples. Every third row is annotated with markings for the peak (or mode) of the distribution (white diamond), as well as the 25th, 50th (median), and 75th percentiles (white tick marks) to highlight the shape of the distribution.</figcaption></figure>

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 (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 ~0.15 dex) are significant relative to the total spread, which is not the case for the [Fe/H] measurements.

### 3.8. Age-Metallicity Relation

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 (AMR) 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 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; Grieves et al. 2018; Lin 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., Feuillet et al. 2019; Hasselquist et al. 2019; Casamiquela et al. 2021; Lian et al. 2022a; 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 AMR 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. 2022a; 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 AMR 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 AMR has less scatter moving toward 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 AMR of the inner galaxy has been more preserved, and less contaminated by migrated stars; or in other words, more stars migrate outward than inward in the Galaxy. This is perhaps not surprising; dynamically, stars in the inner Galaxy migrate outward, and stars in the outer Galaxy migrate inward (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 outward simply because more stars start in the inner Galaxy. The AMR 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 AMR flattens out at young stellar ages, perhaps indicating that chemical
<figure><img src="image_25.png" /><figcaption>**Figure 23.** The age–metallicity relation across the Milky Way disk. Panels represent different spatial zones, laid out in the same way as Figure 8, with rows corresponding to Z and columns increasing in R. The number in the top-right corner of each panel is the number of stars in our sample in that spatial bin. The age and metallicity for individual stars are plotted, colored by [Mg/Fe] abundance. The running median trend is plotted in black square points to guide the eye, with the vertical bars indicating the standard deviation in [Fe/H] for bins in log(age). The typical (median) uncertainty for any given point is shown in the top-right corner of each panel.</figcaption></figure>

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. 2022a). 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 AMR 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 AMR has two separate peaks.

### 3.9. Age–$\alpha$ Relation and Chemical Clocks

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 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 of 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 Figures 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) reported a similar distribution in the AMR 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 (Laporte et al. 2019; Antoja et al. 2020; Ruiz-Lara 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.
<figure><img src="image_26.png" /><figcaption>Figure 24. Same as Figure 23, but showing the contours of point density in each panel.</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>**Figure 25.** The age-$\alpha$ relation across the Milky Way disk. Panels represent different spatial zones, laid out in the same way as Figure 8. The contours represent the density of points on the diagram. The gray background shape outlines the 90% contour for the entire sample, and is the same in all panels for reference. The typical (median) uncertainty for any given point is shown in the top-right corner of each panel.</figcaption></figure>

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 the Appendix to emphasize this point.

### 3.10. Chemical Evolution via Chemical Tagging

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 ≤ [Fe/H] ≤ 0.1 dex). The present-day radial distribution (x-axis) and [Mg/Fe] abundance (y-axis) are 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
<figure><img src="image_28.png" /><figcaption>**Figure 26.** The radial distribution of the low-$\alpha$ (solid lines) and high-$\alpha$ (dashed lines) populations shown as the contour containing 75% of all points for different bins in stellar age (line color) all at fixed $-0.1 \le [Fe/H] \le 0.1$. The center of each contour is marked as a point. The typical (median) uncertainty $[Mg/Fe]$ is plotted in the bottom-left corner.</figcaption></figure>

median value of $[Mg/Fe]$ decreases from $[Mg/Fe] = 0.05$ in the oldest age bin ($10 \le \log(\text{age}) \le 10.1$) to $[Mg/Fe] = 0.00$ in the youngest age bin ($9 \le \log(\text{age}) \le 9.1$). This tracks the chemical evolution of a location in the Galaxy, as Type Ia supernovae begin “diluting” the ISM with iron, thereby decreasing the overall $[Mg/Fe]$ ratio over time.

# 4. Discussion

Using large samples of stars to map the Milky Way in different parameter spaces using metallicity, $\alpha$-element abundances, and age, as demonstrated here, has the potential to place strong constraints on chemical evolution models and reveal the major processes that shaped our Galaxy. Directly comparing our results with specific chemical evolution models is beyond the scope of this paper, but in this discussion section we qualitatively compare our results with predictions from the leading classes of chemical evolution models discussed in Section 1; the “two-infall,” “superposition,” and “clumpy formation” scenarios.

The underlying assumption necessary to interpret these results is that in a well-mixed ISM, stars formed at the same time and the same place in the Galaxy will have the same chemical abundances (both metallicity and $\alpha$-elements). Under this assumption, a spread in abundance at present day for stars at a given age at the same location, whether bimodal or not, can only be produced if stars have moved away from their birth location.

## 4.1. Superposition and Radial Migration

The “superposition” class of evolution models (e.g., Schönrich & Binney 2009a, 2009b; Minchev et al. 2013, 2014, 2017; Johnson et al. 2021) explains 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. 2022a). If stars are formed in situ with a steep metallicity gradient, that gradient will flatten over time as metal-rich inner Galaxy stars migrate outward and metal-poor outer Galaxy stars migrate inward, 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 MDF 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., Roškar et al. 2008; Schönrich & Binney 2009a, 2009b; 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 & Monteiro 2019; Khoperskov et al. 2020b). Using Gaia data, Khoperskov et al. (2020a) estimated the OLR to be located at a Galactocentric radius of around 9 kpc. Khoperskov et al. (2020b) used a high-resolution N-body simulation to investigate the relationship between the OLR and radial migration, and found that stars from the inner Galaxy migrating outward 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 buildup of metal-rich stars at $R \sim 9$ kpc in Figures 5 and 6.

We also see signs of radial migration in the AMR shown in Figure 23. Around the solar neighborhood, there is significant scatter about the age-metallicity trend. If the ISM 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; Minchev et al. 2013; Wang & Zhao 2013; Lian et al. 2022a). While any location in the Galaxy should start with a tight AMR, 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 AMR at larger vertical heights, an effect also seen in our data.
The turnover in the AMR, 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 outward to where they are found today, contaminating the AMR. 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 explored the origin of such a turnover using cosmological simulations, and found that even when the turnover in the AMR 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 the occurrence of stars from the outer Galaxy migrating inward is expected to be a less frequent event 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 midplane 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; Bird et al. 2013; Freudenburg et al. 2017; Bird et al. 2021). Using a quantity like guiding radius ($R_{\text{guide}}$) and maximum height $|Z_{\text{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 (2021), and Spitoni et al. (2022), 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 \le \text{age} \le 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. (2022a) 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 outward instead of inward), 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.

## 4.2. Two-infall or Major Merger

The “two-infall” class of evolution models (e.g., Chiappini et al. 1997, 2001; Spitoni et al. 2019a; Lian et al. 2020a, 2020b, 2020c; Spitoni et al. 2020, 2021, 2022) suggests 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, 2009b; 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) reported 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 outward than inward, 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 reported 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. (2023) found 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 AMR in Figure 23 produces some trends that are better explained by the two-infall model than by radial migration. Minchev et al. (2013, 2014) reported that while the scatter around the AMR 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 AMR 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 AMR 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) showed a disjointed AMR 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 AMR as a function of Galactic position (Figure 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) reported, possibly linked to the most recent interaction with the Sagittarius dwarf spheroidal galaxy as it passed through the disk (e.g., Laporte et al. 2019; Antoja et al. 2020; Ruiz-Lara et al. 2020).

In the age–[Mg/Fe] relation, the age overlap between the high-$\alpha$ and low-$\alpha$ sequences is 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 detected an overlap in age, and claimed 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 the Appendix, 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 multiphase star formation history, such as those presented in the two-infall model, better predicts the trends observed in the inner Galaxy.

### 4.3. Clumpy Star Formation Models

The “clumpy star formation” models (e.g., Clarke et al. 2019; Amarante 2020; Beraldo e Silva 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 Figures 23 and 25, where stars at around log(age) ~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, their Figure 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 have only recently become achievable observationally (e.g., Martig et al. 2021).

## 5. Conclusions

The large sample size, extensive spatial coverage, more complete sampling of the inner Galaxy and midplane inaccessible to optical spectroscopic surveys, and precise abundance measurements for stars in the final data release of APOGEE can help provide strong constraints on Galactic formation and evolution models, in particular for its disk populations. We present results from the final data release of the combined SDSS APOGEE and APOGEE-2 surveys that explore the chemical and age structure of the Milky Way’s disk, measure gradients and distribution functions, and link these new observational constraints to predictions from different chemical evolution models.

Our main conclusions are as follows:
1. Cartography: Overall maps of the Milky Way disk exhibit negative radial age and metallicity gradients. The bar/bulge stands out as more metal-rich and α-poor in the inner Galaxy compared to stars at similar radii but different azimuthal angles (e.g., Wegg et al. 2019; Zasowski et al. 2019; Hasselquist et al. 2020; Eilers et al. 2022).

2. [Mg/Fe] Distribution: The distribution of α-element abundances reveals the chemically bimodal disk structure in the Milky Way (e.g., Fuhrmann 1998; Bensby et al. 2005; Reddy et al. 2006; Lee et al. 2011a; Kordopatis et al. 2015a; Hayden et al. 2015). The low-α disk is thinner (in Z) and more radially extended than the high-α disk (e.g., Yoshii 1982; Gilmore & Reid 1983; Bensby et al. 2011; Bovy et al. 2016). The locus of the low-α sequence varies with radius. Stellar ages within the low-α sequence do not seem to track a true single sequence, implying a superposition of evolutionary tracks.

3. Azimuthal Variance of Metallicity: We find no significant evidence of large-scale azimuthal asymmetry in most of the disk, although the Galactic bar stands out as metal-rich in the mid-height plane. In the solar neighborhood, we see some coherent, nonaxisymmetric structure in metallicity, although it does not obviously correlate with the spiral arms as it does in some studies (e.g., Inno et al. 2019; Poggio et al. 2022; Hawkins 2023).

4. Metallicity Gradients: The Milky Way's full radial metallicity gradient is flat near the center of the Galaxy, and steepens farther out in radius. The high-α disk displays a nearly flat metallicity profile everywhere in the Galaxy, and the low-α disk has a negative gradient that is shallower at high Z than it is close to the plane. We measure the overall radial metallicity gradient of the disk $R \ge 7$ kpc to be $-0.056 \pm 0.001 \text{ dex kpc}^{-1}$. The overall vertical metallicity gradient of the disk at the solar neighborhood is $-0.296 \pm 0.01 \text{ dex kpc}^{-1}$. Both the radial and vertical metallicity gradients flatten with increasing stellar age. These values are consistent with previous studies (e.g., Hartkopf & Yoss 1982; Anders et al. 2014; Hayden et al. 2014; Frankel et al. 2019; Katz 2021; Vickers et al. 2021; Gaia Collaboration et al. 2023).

5. MDF: The MDF of the inner Galaxy has the widest spread, but this narrows with radius. The shape of the MDF skews strongly for the low-α disk, transitioning around $R \sim 9.4$ kpc from having a metal-poor tail in the inner Galaxy to having a metal-rich tail in the outer Galaxy (e.g., Anders et al. 2014; Hayden et al. 2015; Kordopatis et al. 2015b; Loebman et al. 2016; Katz 2021).

6. Age Gradients: Like the metallicity gradient, the age profile of the disk is flat in the inner Galaxy but transitions to a negative gradient in the outer Galaxy. The outer Galaxy's gradient is steeper at higher Z for the low-α population, and flat everywhere for the high-α stars (e.g., Bergemann et al. 2014; Martig et al. 2016; Katz 2021; Anders et al. 2023).

7. ADF: The ADF for the low-α disk changes in skewness similar to the MDF (e.g., Katz 2021), with the inner Galaxy skewed toward younger ages, and the outer Galaxy skewed toward older ages. Above the plane ($|Z| > 1$ kpc), there is significant overlap between the ADF

of the low-α and high-α populations (e.g., Haywood et al. 2013; Hayden et al. 2017; Aguirre et al. 2018; Gent et al. 2022), which does not hold closer to the plane.

8. Age-Metallicity Relation: The AMR exhibits significant spread near the solar neighborhood, but is more tightly constrained in the inner Galaxy and at larger vertical heights (e.g., Casagrande et al. 2011; Bergemann et al. 2014; Feuillet et al. 2018). The slope of the AMR varies with radius, and there exists a population of older, metal-rich stars around the solar neighborhood that are likely present due to radial migration (e.g., Hasselquist et al. 2019; Lian et al. 2022a; Sahlholdt et al. 2022). The AMR suggests that the outer disk began forming low-α stars while the high-α sequence was still forming in the inner disk (e.g., Haywood et al. 2013; Aguirre et al. 2018; Gent et al. 2022).

9. Age-α Relation and Chemical Clocks: The age-α relation appears bimodal nearly everywhere in the Galaxy. The low-α sequence evolves significantly with Galactic position, while the high-α sequence displays a constant trend independent of Galactic position (e.g., Haywood et al. 2013; Feuillet et al. 2018; Katz 2021). There may be evidence of a three-phase star formation history (e.g., Sahlholdt et al. 2022) just outside the solar neighborhood ($9 < R < 12$ kpc).

Our results suggest that radial migration is an important process in shaping the present-day appearance of the disk, especially at large radii and close to the Galactic plane. However, stellar migration alone cannot explain the bimodal nature of the α-element abundances or the distribution of stellar ages in the disk. A noncontinuous evolution model, such as the two-infall scenario or clumpy star formation, appears necessary to explain the trends seen in the inner Galaxy.

## Acknowledgments

We would like to thank the anonymous reviewer for providing useful feedback toward improving this manuscript.

J.I., C.P., J.A.H., and A.S.M. gratefully acknowledge support from NSF grant AST-1909897. D.H.W. and J.A.J. gratefully acknowledge support from NSF grant AST-1909841.

T.C.B. acknowledges partial support for this work from grant PHY 14-30152; Physics Frontier Center/JINA Center for the Evolution of the Elements (JINA-CEE), awarded by the US National Science Foundation, and OISE-1927130: The International Research Network for Nuclear Astrophysics (IReNA), awarded by the US National Science Foundation.

J.G.F.-T. gratefully acknowledges the grant support provided by Proyecto Fondecyt Iniciación No. 11220340, and also from ANID Concurso de Fomento a la Vinculación Internacional para Instituciones de Investigación Regionales (Modalidad corta duración) Proyecto No. FOVI210020, and from the Joint Committee ESO-Government of Chile 2021 (ORP 023/2021), and from Becas Santander Movilidad Internacional Profesores 2022, Banco Santander Chile.

C.A.P. is thankful for funding from the Spanish government through grants AYA2014-56359-P, AYA2017-86389-P, and PID2020-117493GB-100.

Funding for the Sloan Digital Sky Survey IV has been provided by the Alfred P. Sloan Foundation, the U.S. Department of Energy Office of Science, and the Participating Institutions. SDSS acknowledges support and resources from
the Center for High-Performance Computing at the University of Utah. The SDSS website is www.sdss.org.

SDSS is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS Collaboration including the Brazilian Participation Group, the Carnegie Institution for Science, Carnegie Mellon University, Center for Astrophysics | Harvard & Smithsonian (CfA), the Chilean Participation Group, the French Participation Group, Instituto de Astrofísica de Canarias, The Johns Hopkins University, Kayli Institute for the Physics and Mathematics of the Universe (IPMU) / University of Tokyo, the Korean Participation Group, Lawrence Berkeley National Laboratory, Leibniz Institut für Astrophysik Potsdam (AIP), Max-Planck-Institut für Astronomie (MPIA Heidelberg), Max-Planck-Institut für Astrophysik (MPA Garching), Max-Planck-Institut für Extraterrestrische Physik (MPE), National Astronomical Observatories of China, New Mexico State University, New York University, University of Notre Dame, Observatório Nacional / MCTI, The Ohio State University, Pennsylvania State University, Shanghai Astronomical Observatory, United Kingdom Participation Group, Universidad Nacional Autónoma de México, University of Arizona, University of Colorado Boulder, University of Oxford, University of Portsmouth, University of Utah, University of Virginia, University of Washington, University of Wisconsin, Vanderbilt University, and Yale University.

# Appendix
## Stellar Ages and Related Caveats

As introduced in Section 2.4, our age estimates are adopted from the *distmass* value added catalog (Stone-Martinez et al. 2023, submitted), which utilizes a neural network trained on the ASPCAP stellar parameters and APOKASC asteroseismic masses to derive stellar mass for all APOGEE stars. The reported uncertainties, accuracy, and potential caveats related to the *distmass* age estimates are discussed in detail in Stone-Martinez et al. (2023, submitted). Here, we explore and discuss how these uncertainties and our particular sample selection may influence our age-related results. In summary, we find that while the general trends stay the same, both the qualitative and quantitative aspects can differ significantly particularly among the oldest stars. We caution against drawing strong conclusions from stellar ages until more precise ages are available for larger samples of stars. Again, we emphasize that this will only impact our figures that include stellar age.

### A.1. Variation in Spatial Sampling

The *distmass* quality cuts described in Section 2.4 notably result in the exclusion of all stars with $[Fe/H] \le -0.7$ from our sample. In Figure A1, we qualitatively explore where in the Galaxy this selection criterion may influence our results. This shows the fraction of stars $f_{distmass}$ in every spatial bin ($\Delta X = \Delta Y = \Delta Z = 0.5$ kpc), which survive the additional *distmass* quality cut. Close to the solar neighborhood ($6 \le R \le 12$ kpc, $|Z| \le 1$ kpc), the majority of stars pass this additional selection criterion ($f_{distmass} \ge 97\%$). Near the center of the Galaxy ($R \le 5$ kpc, $|Z| \le 1$ kpc), $f_{distmass}$ is generally lower but still exceeds 85%. In the outer Galaxy ($R > 15$ kpc) or beyond the plane ($|Z| \ge 1$ kpc), where most of the stars are expected to be metal-poor, there is a sharp drop in $f_{distmass}$, regularly reaching fractions lower than <80%. This

<figure><img src="image_29.png" /><figcaption>**Figure A1.** The same as Figure 2, but colored by the fraction of stars in our sample that survive the additional distmass quality cuts described in Section 2.4. The purpose of this figure is to qualitatively explore how the distmass criterion might bias our results using stellar age estimates.</figcaption></figure>

shows that our sample with *distmass* ages, while generally a majority, does not equally represent the full sample everywhere in the Galaxy. Particularly beyond $R > 15$ kpc, we caution against directly comparing our results including age to those made with the full sample.

### A.2. Selection of Age Catalog

Neural network-derived ages like *distmass* heavily rely on the accuracy of the training set labels that transfer to the model. Small changes the training set can result in differences in the derived ages of the full training set, which contributes to the large uncertainty values around stellar ages.

The *distmass* catalog publishes results from six different neural network models. All six models were trained using the same subset of stars and the same ASPCAP parameters as training labels, with the only difference being the choice of APOKASC asteroseismic masses used in the training set. The APOKASC catalog publishes asteroseismic results from three different groups (initials “SS,” “MO,” and “TW”—see M.H. Pinsonneault et al. 2023, in preparation for details), referred to as the “uncorrected” masses. Each set of masses also has a “corrected” counterpart, using Gaia parallaxes to calibrate the stellar radius derived from asteroseismology (e.g., Zinn et al. 2019). The calibrations are dependent on stellar log g and are less reliable for stars log g < 2, which is why we elect to use the uncalibrated age set for our final results. A. Stone-Martinez et al. (2023, in preparation) go into more detail on this motivation.
<figure><img src="image_30.png" /><figcaption>**Figure A2.** Systematic uncertainties in the stellar age estimates, comparing the distmass (green and blue) and astroNN (red) catalogs across different ranges of stellar log g (different panels). In each panel, the x-axis is the stellar age from asteroseismology estimates (APOKASC), and the y-axis is the predicted age recovered from the neural network methodology. The typical uncertainty of each point is depicted in black in the upper-left corner of each panel, with the asteroseismic measurements notably having larger uncertainties at lower log g. In both catalogs and across all log g, young stars are generally recovered too old, and old stars are recovered too young. However, there are higher uncertainties in the asteroseismic measurements for the low-log g stars, meaning the neural network age estimates may be more precise than the asteroseismic age estimates in this range.</figcaption></figure>

<figure><img src="image_31.png" /><figcaption>**Figure A3.** The spatial distribution (X-Y face-on view) of stars in APOGEE for different bins of log g (columns). Increasing the log g range of our sample, which would minimize some of the age-related systematics, would also severely limit the area that our sample can probe.</figcaption></figure>

Figure A2 compares different three different age catalogs—the uncalibrated *distmass* ages used in this paper (green points), their calibrated counterparts (blue points), and the *astroNN* APOGEE VAC (Leung & Bovy 2018; Mackereth et al. 2019; red points)—against the APOKASC derived stellar ages on the x-axis for the subset of stars in the *distmass* training set. The comparison is divided by ranges in stellar log g (columns, increasing left to right). We find that all three catalogs suffer a “compression” in ages compared to the APOKASC values, with young stars tending to be assigned ages that are too old, and old stars tending to be fit too young. All three sets of ages predict the ages of higher log g with better accuracy and less scatter. For lower log g stars, the calibrated *distmass* ages predict the young ages slightly better, but *astroNN* and the uncalibrated *distmass* predict the older ages better. For higher log g stars, *astroNN* predicts the young ages slightly better, and *distmass* predicts the older ages better. This offset is thought to be caused by scatter within the training set that the neural network learns from. The asteroseismic masses are more uncertain for low log g stars, which results in a larger mismatch between the APOKASC results and the neural network results. The neural network essentially fits for the relation between [C/N] abundance and stellar age, which should be consistent across a broad range of log g. Therefore, we emphasize for low log g stars that the

neural network likely fits more precise ages than available in the pure asteroseismology results used in the training set.

Although the age estimates show less scatter for higher log g stars, redefining our sample would severely limit the spatial range of the Galaxy we can access. Figure A3 shows the face-on distribution of stars in the disk for various log g limits. Above log g ≥ 2, the outer disk (R ≥ 20 kpc) and the inner disk (R ≤ 5 kpc) disappear from the sample.

Nevertheless, we did test how redefining the sample may influence our results. In general, the observed trends and resulting conclusions do not change, motivating the use of the original sample for the increased spatial range. The radial metallicity gradient as a function of age (Figure 17) differs by less than 0.001 dex kpc⁻¹, and the radial age gradient (Figure 20) differs by 0.03 Gyr kpc⁻¹ using the higher log g sample.

### A.3. Age-α Relation for Different Catalogs

To demonstrate how a different selection of age catalog may influence our results, Figure A4 shows a recreation of Figure 25 for three different age samples, as follows:

1. Test Sample A (log g range): restricting to 2.0 ≤ log g ≤ 2.3, defining a new sample with higher surface gravity;
<figure><img src="image_32.png" /><figcaption>**Figure A4.** Reproduction of the age-α relation (Figure 25) using different samples of stellar ages.</figcaption></figure>
<figure><img src="image_33.png" /><figcaption>Figure A5. Reproduction of the radial age gradients (Figure 20) using different samples of stellar ages.</figcaption></figure>

<figure><img src="image_34.png" /><figcaption>Figure A6. Reproduction of the age-metallicity relation (Figure 23) using different samples of stellar ages.</figcaption></figure>

2. Test Sample B (alternate distmass ages): using the Gaia-corrected APOKASC ages training set from distmass, with the original sample of $1.0 \le \log g \le 2.0$; and

3. Test Sample C (alternate astroNN ages): using astroNN ages instead of the distmass catalog, with the original sample of $1.0 \le \log g \le 2.0$.
Test Sample A (Figure A4 (b)) changes our limits to 2.0 ≤ log g ≤ 2.3, shifting our sample to a different log g range where the age estimates are expected to be more reliable, while purposely avoiding the red clump. However, this suffers from the previously shown lack of spatial distribution, with significantly fewer stars beyond R > 12 kpc and within R < 6 kpc than the original sample used in this paper. This test sample does include a larger number of stars near the solar neighborhood, however. The general location and trends across galactic location are similar, but this sample has more scatter, which leads to larger contour shapes. The bimodality within the low-α sequence (possible evidence for a “third infall”) around 9 < R < 12 kpc is not obvious in this sample, although it is perhaps present at shorter radii (3 < R < 9 kpc).

Test Sample B uses the *distmass* ages using the calibrated ASPCAP training set. The Gaia-motivated calibration is dependent on log g (e.g., Zinn et al. 2019; M.H. Pinsonneault et al. 2023, in preparation), which results in lower log g stars assigned younger ages in the calibrated catalog and higher log g stars being assigned systematically older ages and resulting in an overall age compression in the sample obvious is Figure A4(c). Otherwise, the results are generally similar, with the possible bimodality within the low-α sequence at 9 < R < 12 kpc still present in this sample. The high-α sequence has more spread in age, possibly due to the calibrations propagating uncertainties in the adopted calibrations (from the $T_{eff}$ assumptions, for example) through into the final age estimates.

Test Sample C uses astroNN ages instead of distmass in Figure A4(d), which makes obvious the overall “compression” of the ages caused by the systematics discussed earlier; astroNN ages do not go as old or as young as the distmass ages. The general shapes and trends across the Galaxy otherwise stay the same, and the bimodality within the low-α sequence is not obvious in this sample either.

We also include a reproduction of the median radial age profile (Figure 20) and the AMR (Figure 23) using the higher-log g Test Sample A in Figures A5–A6. As before, the trends in each plot are similar to their low-log g counterparts, but with higher scatter, as the range of Galaxy covered by the sample is more limited despite the lower uncertainties in stellar ages.

### A.4. Summary of Age Caveats

Stellar ages are one of the most powerful quantities available in modern astronomy for revealing the history of the Milky Way, but remain challenging to measure robustly. In this study and others, it is important to consider age-related results with caution and understand the subsequent effects that changes in stellar ages may have on the conclusions. While different age catalogs show generally similar trends for the results of this paper, the quantitative aspects can differ significantly and should be regarded with healthy skepticism until more precise age estimates are available for large samples of stars.

### ORCID iDs

Julie Imig https://orcid.org/0000-0003-2025-3585

Cathryn Price https://orcid.org/0009-0001-4527-6436

Jon A. Holtzman https://orcid.org/0000-0002-9771-9622

Alexander Stone-Martinez https://orcid.org/0000-0003-4761-9305

Steven R. Majewski https://orcid.org/0000-0003-2025-3147

David H. Weinberg https://orcid.org/0000-0001-7775-7261

Jennifer A. Johnson https://orcid.org/0000-0001-7258-1834

Carlos Allende Prieto https://orcid.org/0000-0002-0084-572X

Rachael L. Beaton https://orcid.org/0000-0002-1691-8217

Timothy C. Beers https://orcid.org/0000-0003-4573-6233

Dmitry Bizyaev https://orcid.org/0000-0002-3601-133X

Michael R. Blanton https://orcid.org/0000-0003-1641-6222

Joel R. Brownstein https://orcid.org/0000-0002-8725-1069

Katia Cunha https://orcid.org/0000-0001-6476-0576

José G. Fernández-Trincado https://orcid.org/0000-0003-3526-5052

Diane K. Feuillet https://orcid.org/0000-0002-3101-5921

Sten Hasselquist https://orcid.org/0000-0001-5388-0994

Christian R. Hayes https://orcid.org/0000-0003-2969-2445

Henrik Jönsson https://orcid.org/0000-0002-4912-8609

Richard R. Lane https://orcid.org/0000-0003-1805-0316

Jianhui Lian https://orcid.org/0000-0001-5258-1466

Szabolcs Mészáros https://orcid.org/0000-0001-8237-5209

David L. Nidever https://orcid.org/0000-0002-1793-3689

Annie C. Robin https://orcid.org/0000-0001-8654-9499

Matthew Shetrone https://orcid.org/0000-0003-0509-2656

Verne Smith https://orcid.org/0000-0002-0134-2024

John C. Wilson https://orcid.org/0000-0001-7828-7257

### References

Abdurou'uf, Accetta, K., & Aerts, C. 2022, ApJS, 259, 35

Aguirre, V. S., Bojsen-Hansen, M., Slumstrup, D., et al. 2018, MNRAS, 475, 5487

Allende Prieto, C., Beers, T. C., Wilhelm, R., et al. 2006, ApJ, 636, 804

Amarante, J. A. S., Beraldo e Silva, L., Debattista, V. P., & Smith, M. C. 2020, ApJL, 891, L30

Anders, F., Chiappini, C., Santiago, B. X., et al. 2014, A&A, 564, A115

Anders, F., Chiappini, C., Santiago, B. X., et al. 2018, A&A, 619, A125

Anders, F., Gispert, P., Ratcliffe, B., et al. 2023, arXiv:2304.08276

Antoja, T., Ramos, P., Mateu, C., et al. 2020, A&A, 635, L3

Barbanis, B., & Woltjer, L. 1967, ApJ, 150, 461

Beaton, R. L., Oelkers, R. J., Hayes, C. R., et al. 2021, AJ, 162, 302

Bedell, M., Bean, J. L., Meléndez, J., et al. 2018, ApJ, 865, 68

Bensby, T., Alves-Brito, A., Oey, M. S., Yong, D., & Meléndez, J. 2011, ApJL, 735, L46

Bensby, T., Feltzing, S., Lundström, I., & Ilyin, I. 2005, A&A, 433, 185

Bensby, T., Feltzing, S., & Oey, M. S. 2014, A&A, 562, A71

Bensby, T., Zenn, A. R., Oey, M. S., & Feltzing, S. 2007, ApJL, 663, L13

Beraldo e Silva, L., Debattista, V. P., Khachaturyants, T., et al. 2020, MNRAS, 492, 4716

Bergemann, M., Ruchti, G. R., Serenelli, A., et al. 2014, A&A, 565, A89

Bilitewski, T., & Schönrich, R. 2012, MNRAS, 426, 2266

Bird, J. C., Kazantzidis, S., Weinberg, D. H., et al. 2013, ApJ, 773, 43

Bird, J. C., Loebman, S. R., Weinberg, D. H., et al. 2021, MNRAS, 503, 1815

Bland-Hawthorn, J., & Gerhard, O. 2016, ARA&A, 54, 529

Blanton, M. R., Bershady, M. A., Abolfathi, B., et al. 2017, AJ, 154, 28

Boeche, C., Siebert, A., Piffl, T., et al. 2013, A&A, 559, A59

Bournaud, F., Elmegreen, B. G., & Elmegreen, D. M. 2007, ApJ, 670, 237

Bovy, J., Rix, H.-W., & Hogg, D. W. 2012a, ApJ, 751, 131

Bovy, J., Rix, H.-W., Liu, C., et al. 2012b, ApJ, 753, 148

Bovy, J., Rix, H.-W., Schlafly, E. F., et al. 2016, ApJ, 823, 30

Bowen, I. S., & Vaughan, A. H., Jr. 1973, ApOpt, 12, 1430

Bregman, J. N. 1980, ApJ, 236, 577

Brook, C. B., Kawata, D., Gibson, B. K., & Freeman, K. C. 2004, ApJ, 612, 894

Buder, S., Asplund, M., Duong, L., et al. 2018, MNRAS, 478, 4513

Buder, S., Lind, K., Ness, M. K., et al. 2021, MNRAS, 510, 2407

Carrell, K., Chen, Y., & Zhao, G. 2012, AJ, 144, 185

Casagrande, L., Schönrich, R., Asplund, M., et al. 2011, A&A, 530, A138

Casamiquela, L., Soubiran, C., Jofré, P., et al. 2021, A&A, 652, A25

Chen, B., Hayden, M. R., Sharma, S., et al. 2023, MNRAS, 523, 3791

Cheng, J. Y., Rockosi, C. M., Morrison, H. L., et al. 2012, ApJ, 746, 149
Cheng, X., Anguiano, B., Majewski, S. R., et al. 2020, ApJ, 905, 49

Chiappini, C., Matteucci, F., & Gratton, R. 1997, ApJ, 477, 765

Chiappini, C., Matteucci, F., & Romano, D. 2001, ApJ, 554, 1044

Choi, J., Dotter, A., Conroy, C., et al. 2016, ApJ, 823, 102

Chrobáková, Ž., Nagy, R., & López-Corredoira, M. 2022, A&A, 664, A58

Ciucă, I., Kawata, D., Miglio, A., Davies, G. R., & Grand, R. J. J. 2021, MNRAS, 503, 2814

Clarke, A. J., Debattista, V. P., Nidever, D. L., et al. 2019, MNRAS, 484, 3476

da Silva, R., Porto de Mello, G. F., Milone, A. C., et al. 2012, A&A, 542, A84
Dalcanton, J. J. 2007, ApJ, 658, 941

Dias, W. S., Monteiro, H., Lé pine, J. R. D., & Barros, D. A. 2019, MNRAS, 486, 5726

Donor, J., Frinchaboy, P. M., Cunha, K., et al. 2020, AJ, 159, 199

Edvardsson, B., Andersen, J., Gustafsson, B., et al. 1993, A&A, 275, 101

Eggen, O. J., Lynden-Bell, D., & Sandage, A. R. 1962, ApJ, 136, 748

Eilers, A.-C., Hogg, D. W., Rix, H.-W., et al. 2022, ApJ, 928, 23

Elmegreen, B. G., Elmegreen, D. M., Vollbach, D. R., Foster, E. R., & Ferguson, T. E. 2005, ApJ, 634, 101

Feuillet, D. K., Bovy, J., Holtzman, J., et al. 2018, MNRAS, 477, 2326

Feuillet, D. K., Frankel, N., Lind, K., et al. 2019, MNRAS, 489, 1742

Finlator, K., & Davé, R. 2008, MNRAS, 385, 2181

Frankel, N., Rix, H.-W., Ting, Y.-S., Ness, M., & Hogg, D. W. 2018, ApJ, 865, 96

Frankel, N., Sanders, J., Rix, H.-W., Ting, Y.-S., & Ness, M. 2019, ApJ, 884, 99

Frankel, N., Sanders, J., Ting, Y.-S., & Rix, H.-W. 2020, ApJ, 896, 15

Fraternali, F. 2017, in Gas Accretion onto Galaxies, ed. A. Fox & R. Davé, 430 (Cham: Springer), 323

Freeman, K., & Bland-Hawthorn, J. 2002, ARA&A, 40, 487

Freudenburg, J. K. C., Weinberg, D. H., Hayden, M. R., & Holtzman, J. A. 2017, ApJ, 849, 17

Fuhrmann, K. 1998, A&A, 338, 161

Gaia Collaboration, Brown, A. G. A., & Vallenari, A. 2018, A&A, 616, A1

Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2021, A&A, 650, C3

Gaia Collaboration, Recio-Blanco, A., & Kordopatis, G. 2023, A&A, 674, A38

García Pérez, A. E., Allende Prieto, C., Holtzman, J. A., et al. 2016, AJ, 151, 144

Gebek, A., & Matthee, J. 2022, ApJ, 924, 73

Genovali, K., Lemasse, B., Bono, G., et al. 2014, A&A, 566, A37

Gent, M. R., Eitner, P., Laporte, C. F. P., et al. 2022, arXiv:2206.10949

Gibson, B. K., Pilkington, K., Brook, C. B., Stinson, G. S., & Bailin, J. 2013, A&A, 554, A47

Gilmore, G., Randich, S., Worley, C. C., et al. 2022, A&A, 666, A120

Gilmore, G., & Reid, N. 1983, MNRAS, 202, 1025

Gómez, F. A., Minchev, I., O'Shea, B. W., et al. 2012, MNRAS, 423, 3727

Grieves, N., Ge, J., Thomas, N., et al. 2018, MNRAS, 481, 3244

Gunn, J. E., Siegmund, W. A., Mannery, E. J., et al. 2006, AJ, 131, 2332

Guo, Y., Ferguson, H. C., Bell, E. F., et al. 2015, ApJ, 800, 39

Halle, A., Matteo, P. D., Haywood, M., & Combes, F. 2015, A&A, 578, A58

Hänninen, J., & Flynn, C. 2002, MNRAS, 337, 731

Hartkopf, W. I., & Yoss, K. M. 1982, AJ, 87, 1679

Hasselquist, S., Holtzman, J. A., Shetrone, M., et al. 2019, ApJ, 871, 181

Hasselquist, S., Zasowski, G., Feuillet, D. K., et al. 2020, ApJ, 901, 109

Hawkins, K. 2023, MNRAS, 411, 116

Hawkins, K., Jofré, P., Masseron, T., & Gilmore, G. 2015, MNRAS, 453, 758

Hayden, M. R., Bovy, J., Holtzman, J. A., et al. 2015, ApJ, 808, 132

Hayden, M. R., Holtzman, J. A., Bovy, J., et al. 2014, AJ, 147, 116

Hayden, M. R., Recio-Blanco, A., de Laverny, P., Mikolaitis, S., & Worley, C. C. 2017, A&A, 608, L1

Hayden, M. R., Sharma, S., Bland-Hawthorn, J., et al. 2020, MNRAS, 517, 5325

Haywood, M., Lehnert, M. D., Matteo, P. D., et al. 2016, A&A, 589, A66

Haywood, M., Matteo, P. D., Lehnert, M., et al. 2018, A&A, 618, A78

Haywood, M., Matteo, P. D., Lehnert, M. D., Katz, D., & Gómez, A. 2013, A&A, 560, A109

Haywood, M., Snaith, O., Lehnert, M. D., Matteo, P. D., & Khoperskov, S. 2019, A&A, 625, A105

Helmi, A., Babusiaux, C., Koppelman, H. H., et al. 2018, Natur, 563, 85

Holtzman, J. A., Hasselquist, S., Shetrone, M., et al. 2018, AJ, 156, 125

Holtzman, J. A., Shetrone, M., Johnson, J. A., et al. 2015, AJ, 150, 148

Huang, Y., Schönrich, R., Zhang, H., et al. 2020, ApJS, 249, 29

Iben, I, Jr. 1965, ApJ, 142, 1447

Inno, L., Urbaneja, M. A., Matsunaga, N., et al. 2019, MNRAS, 482, 83

Isern, J. 2019, ApJL, 878, L11

Johnson, J. W., Weinberg, D. H., Vincenzo, F., et al. 2021, MNRAS, 508, 4484

Jönsson, H., Holtzman, J. A., Prieto, C. A., et al. 2020, AJ, 160, 120

Jönsson, H., Prieto, C. A., Holtzman, J. A., et al. 2018, AJ, 156, 126

Jurić, M., Ivezić, Ž., Brooks, A., et al. 2008, ApJ, 673, 864

Katz, D., Gómez, A., Haywood, M., Snaith, O., & Matteo, P. D. 2021, A&A, 655, A111

Kawata, D., & Chiappini, C. 2016, AN, 337, 976

Khoperskov, S., Gerhard, O., Matteo, P. D., et al. 2020a, A&A, 634, L8

Khoperskov, S., Haywood, M., Snaith, O., et al. 2021, MNRAS, 501, 5176

Khoperskov, S., Matteo, P. D., Haywood, M., Gómez, A., & Snaith, O. N. 2020b, A&A, 638, A144

Kobayashi, C., & Nakasato, N. 2011, ApJ, 729, 16

Kordopatis, G., Binney, J., Gilmore, G., et al. 2015b, MNRAS, 447, 3526

Kordopatis, G., Gilmore, G., Wyse, R. F. G., et al. 2013, MNRAS, 436, 3231

Kordopatis, G., Wyse, R. F. G., Gilmore, G., et al. 2015a, A&A, 582, A122

Kubryk, M., Prantzos, N., & Athanassoula, E. 2015, A&A, 580, A126

Lacey, C. G. 1984, MNRAS, 208, 687

Laporte, C. F. P., Minchev, I., Johnston, K. V., & Gómez, F. A. 2019, MNRAS, 485, 3134

Larson, R. B. 1976, MNRAS, 176, 31

Lee, Y. S., Beers, T. C., An, D., et al. 2011a, ApJ, 738, 187

Lee, Y. S., Beers, T. C., Prieto, C. A., et al. 2011b, AJ, 141, 90

Lemasle, B., Hajdu, G., Kovtyukh, V., et al. 2018, A&A, 618, A160

Leung, H. W., & Bovy, J. 2018, MNRAS, 483, 3255

Lian, J., Thomas, D., Maraston, C., et al. 2020a, MNRAS, 494, 2561

Lian, J., Thomas, D., Maraston, C., et al. 2020b, MNRAS, 497, 2371

Lian, J., Zasowski, G., Hasselquist, S., et al. 2020c, MNRAS, 497, 3557

Lian, J., Zasowski, G., Hasselquist, S., et al. 2021, MNRAS, 500, 282

Lian, J., Zasowski, G., Hasselquist, S., et al. 2022a, MNRAS, 511, 5639

Lian, J., Zasowski, G., Mackereth, T., et al. 2022b, MNRAS, 513, 4130

Lin, J., Asplund, M., Ting, Y.-S., et al. 2019, MNRAS, 491, 2043

Lin, J., Dotter, A., Ting, Y.-S., & Asplund, M. 2018, MNRAS, 477, 2966

Linden, S. T., Pryal, M., Hayes, C. R., et al. 2017, ApJ, 842, 49

Loebman, S. R., Debattista, V. P., Nidever, D. L., et al. 2016, ApJL, 818, L6

Lu, Y. L., Ness, M. K., Buck, T., & Carr, C. 2022, MNRAS, 512, 4697

Luo, A.-L., Zhao, Y.-H., Zhao, G., et al. 2015, RAA, 15, 1095

Mackereth, J. T., Bovy, J., Leung, H. W., et al. 2019, MNRAS, 489, 176

Mackereth, J. T., Bovy, J., Schiavon, R. P., et al. 2017, MNRAS, 471, 3057

Mackereth, J. T., Crain, R. A., Schiavon, R. P., et al. 2018, MNRAS, 477, 5072

Magrini, L., Coccato, L., Stanghellini, L., Casasola, V., & Galli, D. 2016, A&A, 588, A91

Magrini, L., Vincenzo, F., Randich, S., et al. 2018, A&A, 618, A102

Majewski, S. R., Schiavon, R. P., Frinchaboy, P. M., et al. 2017, AJ, 154, 94

Marinacci, F., Fraternali, F., Nipoti, C., et al. 2011, MNRAS, 415, 1534

Martig, M., Fouesneau, M., Rix, H.-W., & Edvardsson, B., et al. 2012, AJ, 144, 120

Michtchenko, T. A., Vieira, R. S. S., & Barros, D. A. 2016, A&A, 597, A39
Miglio, A., Chiappini, C., Mackereth, J. T., et al. 2021, A&A, 645, A85

Minchev, I., Anders, F., Recio-Blanco, A., et al. 2018, MNRAS, 481, 1645

Minchev, I., Chiappini, C., & Martig, M. 2013, A&A, 558, A9

Minchev, I., Chiappini, C., & Martig, M. 2014, A&A, 572, A92

Minchev, I., Martig, M., Streich, D., et al. 2015, ApJL, 804, L9

Minchev, I., Steinmetz, M., Chiappini, C., et al. 2017, ApJ, 834, 27

Mollá, M., Díaz, Á., Cavichia, I. O., et al. 2018, MNRAS, 482, 3071

Mor, R., Robin, A. C., Figueras, F., Roca-Fàbrega, S., & Luri, X. 2019, A&A, 624, L1

Myers, N., Donor, J., Spoo, T., et al. 2022, AJ, 164, 85

Ness, M., Hogg, D. W., Rix, H. W., et al. 2016, ApJ, 823, 114

Nidever, D. L., Bovy, J., Bird, J. C., et al. 2014, ApJ, 796, 38

Nidever, D. L., Holtzman, J. A., Allende Prieto, C., et al. 2015, AJ, 150, 173

Nissen, P. E. 2015, A&A, 579, A52

Pilkington, K., Few, C. G., Gibson, B. K., et al. 2012, A&A, 540, A56

Pinsonneault, M. H., Elsworth, Y. P., Tayar, J., et al. 2018, ApJS, 239, 32

Poggio, E., Recio-Blanco, A., Palicio, P. A., et al. 2022, A&A, 666, L4

Price-Jones, N., Bovy, J., Webb, J. J., et al. 2020, MNRAS, 496, 5101

Queiroz, A. B. A., Chiappini, C., Perez-Villegas, A., et al. 2021, A&A, 656, A156

Quinn, P. J., Hernquist, L., & Fullagar, D. P. 1993, ApJ, 403, 74

Rahimi, A., Carrell, K., & Kawata, D. 2014, RAA, 14, 1406

Reddy, B. E., Lambert, D. L., & Prieto, C. A. 2006, MNRAS, 367, 1329

Reid, M. J., Menten, K. M., Brunthaler, A., et al. 2019, ApJ, 885, 131
Robin, A. C., Bienaymé, O., Fernández-Trincado, J. G., & Reylé, C. 2017, A&A, 605, A1

Robin, A. C., Bienaymé, O., Salomon, J. B., et al. 2022, A&A, 667, A98

Romano, D., & Starkenburg, E. 2013, MNRAS, 434, 471

Roškar, R., Debattista, V. P., Quinn, T. R., Stinson, G. S., & Wadsley, J. 2008, ApJL, 684, L79

Ruiz-Lara, T., Gallart, C., Bernard, E. J., & Cassisi, S. 2020, NatAs, 4, 965

Sahlholdt, C. L., Feltzing, S., & Feuillet, D. K. 2022, MNRAS, 510, 4669

Salaris, M., & Cassisi, S. 2005, Evolution of Stars and Stellar Populations (New York: VCH)

Santana, F. A., Beaton, R. L., Covey, K. R., et al. 2021, AJ, 162, 303

Schlesinger, K. J., Johnson, J. A., Rockosi, C. M., et al. 2014, ApJ, 791, 112

Schönrich, R., & Binney, J. 2009a, MNRAS, 396, 203

Schönrich, R., & Binney, J. 2009b, MNRAS, 399, 1145

Sellwood, J. A., & Binney, J. J. 2002, MNRAS, 336, 785

Shapiro, P. R., & Field, G. B. 1976, ApJ, 205, 762

Sharma, S., Hayden, M. R., & Bland-Hawthorn, J. 2021a, MNRAS, 507, 5882

Sharma, S., Hayden, M. R., Bland-Hawthorn, J., et al. 2021b, MNRAS, 510, 734

Shetrone, M., Bizyaev, D., Lawler, J. E., et al. 2015, ApJS, 221, 24

Smith, V. V., Bizyaev, D., Cunha, K., et al. 2021, AJ, 161, 254

Soubiran, C., Bienaymé, O., & Siebert, A. 2003, A&A, 398, 141

Spitoni, E., Aguirre, V. S., Matteucci, F., Calura, F., & Grisoni, V. 2019a, A&A, 623, A60

Spitoni, E., Børsen-Koch, V. A., Verma, K., & Stokholm, A. 2022, A&A, 663, A174

Spitoni, E., Cescutti, G., Minchev, I., et al. 2019b, A&A, 628, A38

Spitoni, E., Recio-Blanco, A., de Laverny, P., et al. 2023, A&A, 670, A109

Spitoni, E., Verma, K., Aguirre, V. S., et al. 2021, A&A, 647, A73

Spitoni, E., Verma, K., Aguirre, V. S., & Calura, F. 2020, A&A, 635, A58

Spitzer, L. J., & Schwarzschild, M. 1951, ApJ, 114, 385

Steinmetz, M., Matijevič, G., Enke, H., et al. 2020, AJ, 160, 82

Stone-Martinez, A., Holtzman, J. A., Imig, J., et al. 2023, AJ, submitted

Thomas, D., Maraston, C., Bender, R., & Mendes de Oliveira, C. 2005, ApJ, 621, 673

Ting, Y.-S., Conroy, C., & Rix, H.-W. 2015, ApJ, 816, 10

Toth, G., & Ostriker, J. P. 1992, ApJ, 389, 5

Toychuchi, D., & Chiba, M. 2018, ApJ, 855, 104

Twarog, B. A. 1980, ApJ, 242, 242

Vázquez, C. V., Magrini, L., Casali, G., et al. 2022, A&A, 660, A135

Vickers, J. J., Shen, J., & Li, Z.-Y. 2021, ApJ, 922, 189

Vincenzo, F., Matteucci, F., Vattakunnel, S., & Lanfranchi, G. A. 2014, MNRAS, 441, 2815

Vincenzo, F., Spitoni, E., Calura, F., et al. 2019, MNRAS Lett., 487, L47

Vincenzo, F., Weinberg, D. H., Miglio, A., Lane, R. R., & Roman-Lopes, A. 2021a, MNRAS, 508, 5903

Vincenzo, F., Weinberg, D. H., Montalbán, J., et al. 2021b, arXiv:2106.03912

Wang, Y., & Zhao, G. 2013, ApJ, 769, 4

Wegg, C., Rojas-Arriagada, A., Schultheis, M., & Gerhard, O. 2019, A&A, 632, A121

Weinberg, D. H., Andrews, B. H., & Freudenburg, J. 2017, ApJ, 837, 183

Weinberg, D. H., Holtzman, J. A., Hasselquist, S., et al. 2019, ApJ, 874, 102

Weinberg, D. H., Holtzman, J. A., Johnson, J. A., et al. 2022, ApJS, 260, 32

Wenger, T. V., Balser, D. S., Anderson, L. D., & Bania, T. M. 2019, ApJ, 887, 114

Wilson, J. C., Hearty, F. R., Skrutskie, M. F., et al. 2019, PASP, 131, 055001

Wu, Y., Xiang, M., Zhao, G., et al. 2019, MNRAS, 484, 5315

Xiang, M., & Rix, H.-W. 2022, Natur, 603, 599

Yanny, B., Rockosi, C., Newberg, H. J., et al. 2009, AJ, 137, 4377

Yoshii, Y. 1982, PASJ, 34, 365

Zamora, O., García-Hernández, D. A., Allende Prieto, C., et al. 2015, AJ, 149, 181

Zasowski, G., Cohen, R. E., Chojnowski, S. D., et al. 2017, AJ, 154, 198

Zasowski, G., Johnson, J. A., Frinchaboy, P. M., et al. 2013, AJ, 146, 81

Zasowski, G., Schultheis, M., Hasselquist, S., et al. 2019, ApJ, 870, 138

Zinn, J. C., Pinsonneault, M. H., Huber, D., et al. 2019, ApJ, 885, 166