Spatial variations in the Milky Way disc metallicity-age relation

TL;DR

This study maps how the relationship between metallicity and age varies across different regions of the Milky Way disc, revealing the effects of radial migration and vertical structure on stellar populations.

Contribution

It provides the first direct observations of spatial variations in the metallicity-age relation in the Milky Way disc using Gaia data and hierarchical modeling.

Findings

Significant variation in [M/H]-age relation with Galactocentric radius and height.

Measured the present-day metallicity gradient as -0.059 dex/kpc.

Confirmed a vertically flared distribution of young stars in the outer disc.

Abstract

Stellar ages are a crucial component to studying the evolution of the Milky Way. Using Gaia DR2 distance estimates, it is now possible to estimate stellar ages for a larger volume of evolved stars through isochrone matching. This work presents [M/H]-age and [$\alpha$/M]-age relations derived for different spatial locations in the Milky Way disc. These relations are derived by hierarchically modelling the star formation history of stars within a given chemical abundance bin. For the first time, we directly observe that significant variation is apparent in the [M/H]-age relation as a function of both Galactocentric radius and distance from the disc mid-plane. The [M/H]-age relations support claims that radial migration has a significant effect in the plane of the disc. Using the [M/H] bin with the youngest mean age at each radial zone in the plane of the disc, the present-day metallicity gradient is measured to be $-0.059 \pm 0.010$ dex kpc$^{-1}$, in agreement with Cepheids and young field stars. We find a vertically flared distribution of young stars in the outer disc, confirming predictions of models and previous observations. The mean age of the [M/H]-[$\alpha$/M] distribution of the solar neighborhood suggests that the high-[M/H] stars are not an evolutionary extension of the low-$\alpha$ sequence. Our observational results are important constraints to Galactic simulations and models of chemical evolution.

Equations

$$\Delta R=\Delta \text{mbox}[M/H]/gradient,$$

Full text

††pubyear: 2019††pagerange: Spatial variations in the Milky Way disc metallicity-age relation–References

1 Introduction

The age-metallicity relation of the Milky Way disc has long been a focus of Galactic evolution studies (e.g. Twarog, 1980; Edvardsson et al., 1993; Casagrande et al., 2011) as any relation found would place tight constraints on models of Galactic chemical evolution (GCE). Simple ‘closed-box’ models of galactic and stellar evolution dictate that over time the mean metallicity of a stellar population will increase as each generation of stars forms out of gas that has been enriched by previous generations. In reality, the absolute timeline of this metallicity enrichment is strongly dependent on the star formation history (SFH) as well as the amount and composition of gas that is injected into or ejected from the system (e.g. Chiappini, Matteucci & Gratton 1997; Dalcanton 2007; Finlator & Davé 2008). In addition, the motions of stars are perturbed in such a way that over time they end up at a different mean distance to the Galactic center, become increasingly eccentric, and/or gain larger vertical oscillations (e.g. Schönrich & Binney, 2009). Therefore an observational charactization of the full Milky Way disc age-metallicity relation provides a concrete result for potentially complex GCE models and cosmological simulations to replicate.

The main difficultly in studying the Galactic age-metallicity relation has been the age determination of FGK field stars, which are the main targets of large spectroscopic surveys. The stellar metallicity and atmospheric parameters can be determined using spectroscopy, but the age cannot be directly measured. Traditionally, the most common method of determining stellar ages, especially for large samples, has been Bayesian isochrone matching based on the method described in Jørgensen & Lindegren (2005). This involves matching the observed properties of a star to models of stellar evolution to infer the most likely age. According to these models, after the first Gyr the observable surface properties of an FGK star change very little during its lifetime on the main sequence but change significantly as it begins to evolve beyond the core hydrogen burning phase. During the turn-off and subgiant phases, the differences in the main observable properties of stars at a given metallicity and different masses (and therefore ages) are large compared to the spectroscopic measurement uncertainties in those parameters. This makes it possible to determine age with an uncertainty of $\sim 1$ Gyr through isochrone matching (e.g. Edvardsson et al., 1993; Casagrande et al., 2011; Bensby et al., 2014, among others). This is unfortunately not the case on the giant branch where stars of different masses and evolutionary stages can have very similar observable properties.

Recent work has found that certain spectral features, such as Balmer lines (Bergemann et al., 2016) or C/N ratios (Masseron & Gilmore, 2015; Martig et al., 2016), can trace the stellar mass, but this is also not a direct measurement of age. While data-driven and neutral network analyses do provide atmospheric parameter, abundance, and age estimates simultaneously from spectra, the age is mainly constrained by the empirically-derived relationship between the stellar mass and the element abundance information in the spectra (e.g. Ness et al., 2016; Mackereth et al., 2019).

Asteroseismology can produce very precise ages for both dwarf and giant stars (see Gai et al., 2011; Chaplin et al., 2014) for regions of the Galaxy observed by Kepler (Borucki et al., 2010; Koch et al., 2010), K2 (Howell et al., 2014), or CoRoT (Baglin et al., 2006a; Baglin et al., 2006b). However, there currently exist many more high-resolution spectroscopic observations of stars outside these fields for which asteroseismology is out of reach. Feuillet et al. (2016, hereafter F16) show that is it possible to determine ages of giant stars to within $\sim 0.18$ dex in log(age) through isochrone matching if precise distance measurements and high-resolution spectroscopy are available. By hierarchically modelling the star formation history (SFH), F16 and Feuillet et al. (2018, hereafter F18) derive [$\alpha$/M]-age and [M/H]-age relations for a sample of solar neighbourhood giants that are in good agreement with results from solar neighbourhood subgiants.

In recent years, the Milky Way disc age-metallicity relation has been examined using observations from large stellar surveys such as the Geneva-Copenhagen Survey (GCS, Casagrande et al., 2011), the Gaia-ESO Survey (GES, Bergemann et al., 2014), the Large sky Area Multi-Object Fiber Spectroscopic Telescope (LAMOST, Xiang et al., 2017), Galactic Archaeology with HERMES (GALAH, Lin et al., 2018; Buder et al., 2019), the Apache Point Observatory Galactic Evolution Experiment (APOGEE, Anders et al., 2017; Silva Aguirre et al., 2018, F18), and the Multi-object APO Radial Velocity Exoplanet Large-area Survey (MARVELS, Grieves et al., 2018). Interestingly, these studies have produced qualitatively similar age-metallicity relations although they have used different types of stars and different age determination methods. Generally, the local age-metallicity relation shows a large spread in metallicity at any given age, with a flat relation for young and intermediate age stars. The metal-poor stars are consistently older, but most studies find some metal-poor stars with intermediate ages and some old stars with solar metallicities.

High-precision studies with smaller number statistics also find a large spread in metallicity at any given age, despite having small age uncertainties (e.g. Bensby et al., 2014; Haywood et al., 2013; Silva Aguirre et al., 2018; Nissen & Gustafsson, 2018). Such studies support the conclusion that the disc age-metallicity relation is intrinsically scattered and the spread is not solely an artefact of observational errors. In contrast, tight age correlations with [$\alpha$/Fe] or other individual elements have been found (e.g. Nissen, 2015; Bedell et al., 2018, F18). The age-metallicity relation is nevertheless interesting as the range of [M/H] values is larger than the range in [X/Fe] ratios. Additionally, it is important to confirm the spread in the age-metallicity relation throughout the Milky Way, which is difficult with smaller samples.

The large range of ages covered by stars with a single metallicity and apparent lack of evolution in the age-metallicity relation has been the major focus of these studies; however, an equally interesting, and perhaps more diagnostic, feature is that the most metal-rich stars have intermediate ages. This is found in most studies of the solar neighbourhood, but is most striking in Figure 3 of F18, who examine the age distribution of mono-metallicity bins. The spread in age is large in each metallicity bin, but the mean age is nevertheless well defined with this technique and chemical evolution can be clearly seen from old metal-poor stars to younger solar metallicity stars. This trend is reversed at metallicities above solar, producing a turnover feature. The metal-rich stars are on average older than the solar metallicity stars. The hierarchical modelling technique of F16 and F18 relies on determining the mean age of a group of stars (in this case grouped by metallicity), therefore the trends derived are in fact metallicity-age relations (MAR).

Minchev et al. (2013) use a GCE model in a full cosmological context, which includes radial migration, that reproduces these features of the observed local age-metallicity relations. Radial migration is the migration of stars inward or outward from their birth radius while maintaining circular, or near circular, orbits (see Wielen et al., 1996; Sellwood & Binney, 2002; Schönrich & Binney, 2009; Loebman et al., 2011). This stellar migration is caused by the conservation of angular momentum during interactions with resonant features in the disc such as spiral arms. In the model of Minchev et al. (2013), the stars currently at the solar Galactocentric radius were born at a range of radii, see their Figure 3. The metal-rich stars currently in the solar neighborhood were preferentially born in the inner Galaxy where star formation began earlier and proceeded more rapidly, resulting in metal-rich stars born earlier than was possible at larger Galactocentric radii. GCE models that include radial migration can explain the large spread in metallicity at all ages because the birth age-metallicity relation is different at each Galactocentric radius, so as stars migrate inward or outward, they pollute the local age-metallicity relation of their new Galactocentric orbit (e.g. Schönrich & Binney, 2009; Minchev et al., 2013; Kubryk et al., 2015).

Mackereth et al. (2017) find a similar result observationally using APOGEE red giant branch stars and ages based on the mass-CN relation (Martig et al., 2016). They show that mono-age, mono-[Fe/H] stars in the low-$\alpha$  sequence have donut-like surface-mass density profiles. Stars younger than 3 Gyr are more tightly concentrated around the peak density, while the distribution of older stars is broader because they have radially migrated from their birth radius. The average metallicity of the young stars is a function of Galactocentric distance, with the metal-poor stars being concentrated in the outer Galaxy and the metal-rich stars being concentrated in the inner Galaxy. The radial metallicity gradient of the disc has been found by many previous studies using large spectroscopic surveys such as GCS (Casagrande et al., 2011), SEGUE (Lee et al., 2011), RAVE (Boeche et al., 2013, 2014), GES (Bergemann et al., 2014), and APOGEE (Hayden et al., 2014; Anders et al., 2014, 2017). Studies focusing on young stars to measure the gradient have used Cepheids (e.g. Genovali et al., 2014; Inno et al., 2019) and open clusters (e.g. Reddy et al., 2016; Donor et al., 2018). The radial metallicity gradient supports the idea that the metal-rich stars with intermediate ages were likely born in the inner Galaxy.

While radial migration is predicted to have a significant effect in the plane of the disc, the process is less efficient for stars with larger vertical velocities, therefore it should have a smaller effect on stars at larger distances from the plane of the disc. This has been shown in dynamical models (e.g. Solway et al., 2012), and simulations of galactic discs (e.g. Bird et al., 2013). Observationally, Hayden et al. (2015) find that the shape of the metallicity distribution functions for stars through the disc are consistent with a metallicity gradient and a simple model of radial migration in the plane. Moving away from the plane, the metallicity distribution functions are uniform at all Galactocentric radii, perhaps reflecting a homogeneously mixed gas disc rather than a diffusive radial migration process.

Until recently the samples of stars with age determinations have mostly been limited to the local disc or stars with chemically-based ages. But with the extensive parallax measurements provided by the second Data Release of the European Space Agency Gaia Mission (DR2, Gaia Collaboration et al., 2018), it is now possible to derive precise age trends for stars beyond the solar neighborhood. This work explores the spatial variations in the MAR of the Milky Way disc. Using these variations we discuss the potential dependence of radial migration efficiency on scale height, the disc radial metallicity gradient, and the age evolution of the high-$\alpha$  and low-$\alpha$  sequences.

2 Sample and Methods

The sample presented is taken from the Data Release 14 (DR14, Abolfathi et al., 2018) of the Sloan Digital Sky Survey-IV (SDSS-IV Blanton et al., 2017) APOGEE (Majewski et al., 2017) for which Gaia Data Release 2 (GDR2, Gaia Collaboration et al., 2018) parallax measurements are available. APOGEE is an $H$-band spectroscopic survey of Milky Way stars using the 300-fibre fed, high-resolution ($R\sim$ 22,500), APOGEE spectrograph (Wilson et al., 2019) on the Sloan Foundation 2.5-m telescope (Gunn et al., 2006). APOGEE data are reduced using a standard pipeline Nidever et al. (2015) and analyzed by the APOGEE Stellar Parameters and Chemical Abundances Pipeline (ASPCAP, García Pérez et al., 2016). In DR14, ASPCAP uses a specially computed library of 1D plane parallel stellar atmospheric models assuming local thermodynamic equilibrium (Zamora et al., 2015), which is calculated using a custom built line-list (Shetrone et al., 2015), to simultaneously determine effective temperature ($T_{\mbox{\footnotesize eff}}$), surface gravity ($\log g$), metallicity ([M/H]), [$\alpha$/M], [C/M], and [N/M] through $\chi^{2}$-minimization. Calibration and quality analysis of DR14 are described in Holtzman et al. (2018) and Jönsson et al. (2018). The targeting scheme for APOGEE is described in Zasowski et al. (2013) and Zasowski et al. (2017).

The following selection and quality criteria have been applied in this work:

• •

$\sigma_{\pi}/\pi<0.2$

• •

$(J-K)_{0}\geq 0.5$

• •

1.0 < $\log g$< 3.8

• •

3500 < $T_{\mbox{\footnotesize eff}}$< 5500

• •

SNR > 80.

• •

No TAR_BAD or \verb VINI_WARN ASPCAP flag set

• •

No AD_PIXELS or \verb VERY_RIGHT_NEIGHBOR star flag set

• •

No known or candidate cluster members

The final sample contains 81,400 stars. The results presented have removed known cluster members and are limited to stars with radial Galactocentric distances ($R_{\rm Gal}$) between 5 and 13 kpc and a distance from the plane ($|z|$) of less than 2 kpc, resulting in 77,562 stars. To explore the spatial variation in the derived age trends, the sample is divided into four $R_{\rm Gal}$  bins (5-7, 7-9, 9-11, and 11-13 kpc) and three $|z|$  bins (0.0-0.5, 0.5-1.0, and 1.0-2.0 kpc), resulting in 12 disc zones. With increasing $|z|$, the zones become dominated by stars with positive $z$ height, but we find the MARs are symmetric in $z$. We therefore use $|z|$  to increase the signal in all zones. Figure 1 shows the $\log g$-$T_{\mbox{\footnotesize eff}}$  diagram for each of the 12 zones with the [M/H] indicated by the color. The full red giant branch is sampled in all zones except $R_{\rm Gal}$  11-13 kpc where the high $\log g$  stars are too faint for the Gaia parallax restriction.

2.1 Age determination

The age trends presented here were determined using the hierarchical modelling method described in F16 and F18. This method constrains the parameters of a model SFH by combining the age likelihood functions produced from Bayesian isochrone matching for a group of stars. The parameters used in the Bayesian isochrone matching are $T_{\mbox{\footnotesize eff}}$, $\log g$, overall metallicity, and absolute $K$ magnitude. The overall metallicity is calculated from the APOGEE [M/H] and [$\alpha$/M] using the prescription of Salaris et al. (1993). The stellar atmospheric parameters ($T_{\mbox{\footnotesize eff}}$, $\log g$, [M/H], and [$\alpha$/M]) are calibrated APOGEE DR14 values. A full description of the calibrations applied to DR14 is available in Holtzman et al. (2018) and Jönsson et al. (2018) provides a thorough comparison of APOGEE abundances to independent analyses. The absolute $K$ magnitudes are calculated using the $K$-band magnitudes from the Two Micron All Sky Survey (Skrutskie et al., 2006), a distance, and an extinction. The distance value is taken from Bailer-Jones et al. (2018). The extinction is taken from the APOGEE targeting information using the K_TRG   parameter (see Zasowski et al., 2013) which uses the RJCE method (Majewski et al., 2011). If the K_TRG   parameter is not available, then the WISE $K$-band extinction is used from the K_WISE , parameter. In this work, PRSEC isochrones (Bressan et al., 2012) are used with a lognormal Chabrier initial mass function (Chabrier, 2001). We note that PARSEC isochrones do not include atomic diffusion, which can affect the surface metallicity, but this process is thought to have a minimal effect on giant stars (Dotter et al., 2017).

The model SFH used is a Gaussian function plus a uniform component allowing for outliers within a group of stars assuming an outlier fraction of 7.5 % as in F16. For this work, the stars are grouped by abundance, specifically [M/H] and [$\alpha$/M]. The widths of the abundance bins are determined by the mean uncertainty in the abundance measurements. If fewer than 15 stars lie within the bin then the width is increased until it contains 15 stars. The result of the hierarchical modelling is a mean age and age dispersion for the stars in a given abundance bin and can be applied to individual abundances independently. Due to the method of binning stars in abundance starting with the stars with the lowest abundances, the lowest abundance bin usually contains stars that cover a larger abundance range and are typically outliers from the main abundance distribution.

We apply a correction for the bias in $\log g$  imposed by the 20 % cut in parallax uncertainty. The low $\log g$  stars are intrinsically bright and are therefore preferentially farther away in a survey selecting targets from a limited magnitude range. At a given apparent magnitude, more distant stars have larger parallax uncertainties than closer stars. These effects lead to a bias against low $\log g$  stars in this sample. The $\log g$  bias in this sample is illustrated in Figure 2. The top panel shows the $\log g$  distribution of the sample with no parallax uncertainty cut in black, and imposing a maximum of 20 % in parallax uncertainty in red. The bottom panel shows the fraction of stars recovered after imposing the 20 % cut as a function of $\log g$. While the parallax uncertainty limit has almost no effect on the high $\log g$  stars, up to 80 % of the low $\log g$  stars are removed by this selection. This selection bias is accounted for during the normalization of the SFH model, see $N(a)$ in section 4.6 of F16.

A selection function that accounts for the APOGEE-1 color selection and the cuts we imposed on $T_{\mbox{\footnotesize eff}}$  and $\log g$  is also applied during the hierarchical modelling. This does not account for the full APOGEE - Gaia crossmatch selection function, which is more complex and beyond the scope of this paper. A full treatment of the survey selection functions is planned for future comparisons with Galactic models where accounting for possible biases will be crucial. Here we discuss generally the biases expected due to the APOGEE-1 and APOGEE-2 targeting strategies (Zasowski et al., 2013, 2017). APOGEE-1 used a single color selection for the disc and bulge fields, $(J-K)_{0}\geq 0.5$, but APOGEE-2 applies a dual color selection, $0.5\leq(J-K)_{0}\leq 0.8$ and $(J-K)_{0}\geq 0.8$. In both APOGEE-1 and APOGEE-2 the halo fields use a color selection of $(J-K)_{0}\geq 0.3$, but our full sample color cut in this work is also applied to the halo fields, so no age bias is expected. Possible age biases caused by the APOGEE-2 dual color selection are expected to be up to $\sim 0.1$ dex older in the disc and bulge fields.

The most distant zones are biased towards luminous, upper giants branch stars, as seen in Figure 1, due to the parallax uncertainty cut. This luminosity bias is expected to cause a bias towards younger ages, up to $\sim 0.15$ dex, in the most distant zones. Using the solar neighborhood sample, we find that there is no bias in [M/H] or [$\alpha$/M] due to the lack of lower giant branch stars. We therefore determine that the shape of the outer zone abundance-age relations should be unaffected. This luminosity age bias is opposite to the effects expected from the color bias in the disc and bulge fields. In Section 3 we discuss how these biases may affect our interpretation of results.

3 Age trends

3.1 Metallicity-Age Relation

The age-metallicity relation in the solar neighborhood has been observed to have stars with a large spread in metallicity at any given age (e.g. Edvardsson et al., 1993; Casagrande et al., 2011; Bensby et al., 2014). Most recent studies find that the most metal-poor stars have the oldest ages, the most metal-rich stars have intermediate ages, and the youngest stars have solar metallicities (see Casagrande et al., 2011, F18). These observed deviations from the narrow and monotonic age-metallicity relation predicted by simple models of chemical evolution have been tentatively attributed to radial migration of stars in the Galactic disc (F18). With the sample presented here it is possible to search for detailed spatial variations in the MAR through the Galactic disc. Again, we note that the analysis presented determines a mean age for stars binned in [M/H]. Comparisons are only made to unbinned literature age-metallicity relation results.

Figure 3 shows the hierarchically modelled MAR for 12 disc zones. The points mark the mean age derived for each abundance bin and the error bars indicate the uncertainty in the mean age. The blue shaded regions show the age dispersion and abundance width of each bin. Bins containing only 15 stars are given a lighter color. From this figure is it clear that significant spatial variations exist in the Galactic MAR. To put these MARs in context with the [$\alpha$/M] vs [M/H] abundance distributions, Figure 4 shows the [$\alpha$/M] vs [M/H] distribution for each of the 12 zones. These distributions show that the high-$\alpha$  sequence is dominant at high $|z|$  and has a shorter scale length than the low-$\alpha$  sequence, which dominates in the plane of the disc and is strongly present out to larger $R_{\rm Gal}$. The relative spatial distributions of the high- and low-$\alpha$  sequences shown here are in agreement with APOGEE DR12 results from Hayden et al. (2015) and the APOGEE DR14 [Fe/Mg] vs [Mg/H] distributions shown in Weinberg et al. (2019).

Starting in the solar neighborhood, the clear turnover behavior of the local MAR noted by F18 is recovered in this larger sample. This sample contains 15,180 stars in the $7<$ R_{\rm Gal}$$<9, $0<$|z|$<0.5$ zone, while F18 presented 721 stars within 400 pc of the Sun. The increase in mean age with decreasing metallicity below [M/H] $\sim 0$ extends to [M/H] $\sim-1$. The youngest stars in this zone have a metallicity of -0.1 to 0.2, consistent with the local interstellar medium (Nieva & Przybilla, 2012), and most were likely born within the zone. The Sun is older than mean age of solar [M/H] stars in the solar neighborhood. Assuming an age of $4.66$ Gyr (Dziembowski et al., 1999) and [M/H] of 0, Figure 3 suggests that the Sun was likely born inward of its current $R_{\rm Gal}$  between 5 and 7 kpc, or even closer to the Galactic centre. This is consistent with previous estimates (Wielen et al., 1996; Minchev et al., 2013, 2018). The metal-rich stars, which are older than the solar metallicity stars and more evolved chemically, probably radially migrated from elsewhere in the disc. The presence of intermediate-age, metal-rich stars in the solar neighborhood is a generic result of GCE models that include radial migration. In such models, these metal-rich stars form in the inner disc, see figure 5 of Minchev et al. (2014). In Kubryk et al. (2015), this results in a MAR with a similar turnover feature as the present work when considering the mean age of all stars, not just those formed in situ, see their figure 12.

One of the most striking features of Figure 3 is the presence of the mean age turnover of high metallicity stars in all of the $0<$|z|$<0.5$ zones (bottom row of Figure 3). This suggests that the radial migration of stars in the disc could be a significant process (e.g. Minchev et al., 2014; Frankel et al., 2018; Weinberg et al., 2019). In this larger sample, there is also a secondary turnover feature around [M/H] $\sim 0.4$ dex. The cause of this feature is unknown, but could reflect the intrinsic age-metallicity relation of the inner disc if these highest metallicities are significantly dominated by stars migrated from a similar birth $R_{\rm Gal}$. The presence of a second turnover feature is suggested by figure 15 of F18, which uses a mixture of analytic chemical evolution models representing multiple zones of chemical evolution to simulate a population of stars born over a range of $R_{\rm Gal}$. This second turnover is also seen in figure 12 of Kubryk et al. (2015) and is caused by stars above [Fe/H] = 0.4 coming from the inner $2-3$ kpc.

Interestingly, the metallicity of the primary turnover changes as a function of $R_{\rm Gal}$. The metallicity of the youngest stars in each zone is shown by the left-facing triangles in Figure 3. This is a confirmation of the disc radial metallicity gradient predicted by simulations and observed using Cepheids (e.g. Genovali et al., 2014; Inno et al., 2019), open clusters (e.g. Reddy et al., 2016; Donor et al., 2018), and field stars (e.g. Boeche et al., 2013; Bergemann et al., 2014; Hayden et al., 2014). The confirmation of the metallicity gradient supports the hypothesis that the most metal-rich stars in the plane of the disc likely came from the inner Galaxy, R_{\rm Gal}$$<5 kpc, and the assumption used in many GCE models that star formation began earlier and proceeded more rapidly in the inner Galaxy.

Using the youngest mean age and mean $R_{\rm Gal}$  of these four radial zones in the plane of the disc, we estimate the metallicity gradient to be $-0.061\pm 0.015$ dex kpc*-1*. If we use $R_{\rm Gal}$  bins of 1 kpc instead of 2 kpc, the measured metallicity gradient is $-0.059\pm 0.010$ dex kpc*-1*. These two measurements are consistent, but higher precision is reached with a finer $R_{\rm Gal}$  binning. Our gradient measurement is consistent with measurements from Cepheids by Genovali et al. (2014) and Lemasle et al. (2007), recent open cluster measurements by Donor et al. (2018), APOGEE measurements of young field giants with ages determined using [C/N] abundances (Hasselquist et al., 2018) and asteroseismology (Anders et al., 2014, 2017), young field dwarfs observed by GES (Bergemann et al., 2014) and RAVE (Boeche et al., 2013), and the present-day gradient reported by Minchev et al. (2018) using young field subgiant and turn-off stars. However, previous studies using open clusters and giant stars of all ages measure $\sim 0.08-0.1$ dex kpc*-1* (e.g. Frinchaboy et al., 2013; Hayden et al., 2014; Jacobson et al., 2016). This difference could come from increased sample sizes, improved distance measurements, and an emphasis on using young stars in the other studies.

In the $7<$ R_{\rm Gal}$$<9 kpc bin, at larger distances from the disc plane, the turnover softens and almost flattens in the $1<$|z|$<2$ kpc bin where the mean age does not vary much in different metallicity bins. The softening in the turnover with distance from the disc plane could suggest that stellar migration is much more efficient in the plane of the disc and has only a small effect at large $|z|$. As seen in Figure 4, there are also fewer metal-rich stars in the highest $|z|$  zone, again consistent with the picture that the higher $|z|$  zones do not contain many radially migrated stars. However, the change in the turnover feature with increasing $|z|$, could also be caused by a lack of young stars at high $|z|$. The $7<$ R_{\rm Gal}$$<9, $1<$|z|$<2$ bin is dominated by the high-$\alpha$  sequence, which has been found to be uniformly old (see Figure 7 and Xiang et al., 2017).

With increasing $|z|$, the youngest stars are no longer present in the inner disc, R_{\rm Gal}$$<9, and the zones are dominated by older, high-$\alpha$  stars. At these $R_{\rm Gal}$  the young stars are born kinematically cold and close to the plane of the disc. With time their vertical velocities may increase, allowing them to spend time at larger $|z|$  where the older, high-$\alpha$  sequence dominates. The $5<$ R_{\rm Gal}$$<7, $1<$|z|$<2$ zone is almost entirely dominated by the high-$\alpha$  sequence in this sample, and shows very little change in mean age as a function of [M/H]. This suggests that the high-$\alpha$  sequence is either well mixed or formed all at once.

In the disc plane, the mean age of the metallicity bin with the youngest stars decreases as a function of $R_{\rm Gal}$  from 2-3 Gyr in the inner two zones to 1-2 Gyr in the outer two zones. There is an increased dominance of young stars in the outer disc at all $|z|$  heights. This dominance of young stars in the outer Galaxy has been observed using APOGEE (Ness et al., 2016) and LAMOST (Xiang et al., 2017). It is likely related to the decreased relative contribution in the outer disc of the high-$\alpha$  sequence, which has been found to be generally older than the low-$\alpha$  sequence. While we expect some selection effects are present, discussed in Section 2, we estimate that the net age bias is quite small. The low $|z|$  zones could be biased towards older ages, which would only amplify the trends observed with $|z|$. The $11<$ R_{\rm Gal}$$<13 zones would be most significantly affected by the luminosity bias towards younger ages, but we note that the presence of young stars at larger $|z|$  is clearly seen in the $9<$ R_{\rm Gal}$$<11 zones as well.

Overall, the spatial variations in the MAR are in excellent agreement with the age trends found by Hasselquist et al. (2018) inferred using APOGEE [C/N] abundance ratios, see their figure 4. It has been shown that the [C/N] of giant stars correlates with the mass, and therefore the age, due to internal mixing of CN-cycle processed material from the core (Masseron & Gilmore, 2015; Martig et al., 2016). This relation was confirmed using APOGEE DR14 data by Hasselquist et al. (2018) and F18. Hasselquist et al. (2018) find a turnover in the [C/N] at high [Fe/H], corresponding to a turnover in age, in the plane of the disc that weakens with $|z|$. They also find evidence of fewer old stars in the outer disc.

3.2 Flared young disc

In the present sample, there is an increased fraction of young stars at larger distances from the mid-plane in the outer Galaxy compared to the inner Galaxy. In Figure 3, the outer disc shows young stars present at larger $|z|$  zones than in the inner disc, reflecting the flared age distribution of the outer disc observed using stellar ages by Ness et al. (2016) and Xiang et al. (2017), and predicted by simulation (e.g. Minchev et al., 2014; Minchev et al., 2015; Rahimi et al., 2014). Bensby et al. (2011) also note an absence of $\alpha$-rich stars in the outer disc, consistent with both Figure 4 and an increased fraction of young stars in the outer disc.

To better illustrate the flared distribution of the young stars, we hierarchically model the mean age of the sample as a function of $R_{\rm Gal}$  and $|z|$. Figure 5 shows the $|z|$  vs $R_{\rm Gal}$  distribution of the sample colored by the mean age of each bin in the left panel and by the number of stars in each bin in the right panel. The distribution is binned by 200 pc in $R_{\rm Gal}$  and 100 pc in $|z|$. As in Figure 3, bins are required to contain at least 15 stars. In this case the bin size is not increased in order to contain 15 stars, the bin is simply not shown. The mean age of the outer disc is younger than the inner disc. In particular, inwards of 9 kpc, the young stars dominate at $|z|$  less than 300 pc. Beyond $R_{\rm Gal}$  of 9 kpc, the young stars dominate at larger $|z|$, reaching 1 kpc at $R_{\rm Gal}$  of 12 kpc. This is consistent with inside-out disc formation (e.g. Bird et al., 2013), but inspection of the density profiles of mono-age populations is needed to confirm the true spatial distribution of young stars.

The luminosity bias in the outer disc is expected to be counteracted by the color bias in the disc fields. At higher $|z|$  the density of stars is low and not strongly represented in this figure. The flared behavior is also clear even at 9-10 kpc. The flared age structure is suspiciously centered near the solar $R_{\rm Gal}$  (8 kpc) and could be influenced by some selection effects. However, the overall shape is similar to Xiang et al. (2017, see their figure 23) using LAMOST data with different selection criteria. We preformed the same analysis on a sample of only APOGEE-1 stars (for which the single color selection has been taken into account) and the resulting age distribution is the same as Figure 5.

3.3 [$\alpha$/M]-Age Relation

Figure 6 is the same as Figure 3, but binned in [$\alpha$/M] abundance instead on [M/H]. The solar neighborhood [$\alpha$/M]-age relation (AAR) in this sample is in good agreement with the AAR presented in F18. The mean age increases rapidly with increasing [$\alpha$/M] at low abundances, but at [$\alpha$/M] above 0.15 dex, the mean age is approximately constant. Other $R_{\rm Gal}$  zones in the disc plane have a smoother transition in the AAR from low-$\alpha$  to high-$\alpha$  stars bins, but most zones show a transition around [$\alpha$/M] $\sim 0.1$ dex, above which the high-$\alpha$  sequence dominates. This transition in AAR suggests that the high-$\alpha$  sequence and the low-$\alpha$  sequence had very different chemical enrichment histories. The large age evolution at low [$\alpha$/M] and smaller age evolution at high [$\alpha$/M] is consistent with most models of GCE (e.g. Kubryk et al., 2015) and observations of the local Galaxy by Silva Aguirre et al. (2018) using asteroseismology in the Kepler field as well as Haywood et al. (2013) using solar neighbourhood subgiants.

The turnover towards older ages in the lowest abundance bins of the $7<$ R_{\rm Gal}$$<9, $0<$|z|$<0.5$ zone is reminiscent of the age trends with Si, S, and Ca noted by F18. As in F18, the stars populating these abundance bins are mainly outliers in the [$\alpha$/M] vs [M/H] distribution and are not strongly present in other spatial zones. It is possible that these stars have been accreted from a merger, such as Sagittarius, and do not belong to the disc population. The highest [$\alpha$/M] bins in this zone have mean ages that are younger than the bin at 0.2 dex. The cause of this is unknown, but these bins contain fewer stars and have large mean age uncertainties.

The AAR is very similar in all the spatial zones; the main differences arise due to the presence of young stars. This is apparent in the $7<$ R_{\rm Gal}$$<9 zone. At farther distances from the midplane the young stars are no longer present, as noted in the MAR. This results in older mean ages for the low [$\alpha$/M] bins, making the full relation appear steeper. The AAR for the [$\alpha$/M] bins above $\sim 0.1$ is similar in all zones suggesting that the high-$\alpha$  sequence is fairly uniform across the disc. In the outer disc, the young stars are present at larger $|z|$  than in the inner disc, again lending evidence to the flared distribution of young stars in the outer disc. The outer disc reaches mean ages of only 4-5 Gyr. As discussed above, the high-$\alpha$  sequence is not strongly present in the outer disc and it is unlikely that star formation rates were high at these $R_{\rm Gal}$  in the early Milky Way. As in Figure 3, the outer disc zones are most likely to be biased by selection effects, but the $9<$ R_{\rm Gal}$$<11 zones already have an increased presence of young stars.

3.4 [M/H]-[$\alpha$/M]-Age Distribution

Figure 7 combines the age trends presented in Figures 3 and 6 for the solar neighborhood zone, $7<$ R_{\rm Gal}$$<9, $0<$|z|$<0.5$. In the left panel the mean age is represented by the color of the bins in the [$\alpha$/M] vs [M/H] distribution. The right panel is colored by the logarithmic number density of stars. This distribution is binned by 0.05 dex in [M/H] and 0.03 dex in [$\alpha$/M]. Bins with fewer than 15 stars are not shown. The typical uncertainty in the mean age in this figure is 0.09 dex, or 0.2 Gyr at 1 Gyr and 1.7 Gyr at 7 Gyr. Here we provide a qualitative interpretation of this figure in the context of the literature and GCE.

In this figure the youngest stars are concentrated around [M/H] $\sim 0.1$ and [$\alpha$/M] $\sim-0.05$, in agreement with the single-element age relations in Figures 3 and 6. Note that the lowest [$\alpha$/M] stars from Figure 6 that had intermediate ages are not present in high enough numbers to appear in Figure 7. Here it is obvious that the high metallicity stars ([M/H] $>0.3$) have both older ages and higher [$\alpha$/M] than the youngest stars.

The high-$\alpha$  sequence is older and shows little age evolution with [M/H] and [$\alpha$/M] along the sequence until solar metallicity, as is seen in Figure 6. This is consistent with the $\sim 1$ Gyr time delay of type Ia supernovae compared the type II, which causes the ‘knee’ in the high-$\alpha$  sequence (Tinsley, 1979; Matteucci & Greggio, 1986). The high metallicity of the knee, relative to local dwarf galaxies, indicates that the star formation rate was high at that time (see Tolstoy et al., 2009, and references therein). The steeper slope of the high-$\alpha$  sequence compared to the low-$\alpha$  sequence is likely caused by a decline in star formation, which results in a dominant contribution of the delayed type Ia supernovae over type II (see Hill et al., 2019). This is consistent with derivations of the local star formation history that find a burst of significant star formation at early time followed by a lull in star formation using thick disc stars (Snaith et al., 2014; Haywood et al., 2015) and a two-infall model (Spitoni et al., 2019). Most GCE models find the Milky Way disc star formation history can be approximated by a peak in star formation at around 9-10 Gyr ago, followed by a exponential decline in star formation (e.g. Kubryk et al., 2015; Rybizki et al., 2017; Côté et al., 2017), also consistent with our results.

The high [M/H] stars appear to be an extension of the high-$\alpha$  sequence rather than part of the low-$\alpha$  sequence, as debated in Nidever et al. (2014). Again, it is unlikely that these stars formed locally due to their intermediate age and high [M/H]. In Anders et al. (2018), these stars are referred to as Inner Disk III and Inner Disk IV based on significant separation of the stars in a t-distributed stochastic neighbour embedding (t-SNE) analysis and their cold orbits. The low-$\alpha$  sequence is younger and shows more age evolution with [M/H] and [$\alpha$/M] along the sequence than the high-$\alpha$  sequence, confirming the more extended star formation history of the low-$\alpha$  sequence proposed previously (e.g. van Dokkum et al., 2013; Snaith et al., 2014; Rybizki et al., 2017).

The metal-poor end of the low-$\alpha$  sequence is approximately the same age as the metal-rich extension of the high-$\alpha$  sequence. If the high-[M/H] stars formed in the inner disc (Anders et al., 2018) and the low-$\alpha$  sequence formed throughout the disc after some gas accretion event (as suggested by e.g. Chiappini et al., 1997, 2001), then it is likely that the metal-poor, low-$\alpha$  stars formed in the outer disc and the inner disc was not significantly diluted with the accreted material. If the accreted gas did indeed reach the inner disc and the metal-rich stars formed post accretion (coeval with the metal-poor, low-$\alpha$  stars), then inner disc metallicities must have reached [M/H] $>+0.5$ before the gas accretion. Very few stars are have been reported in the literature with such high metallicities (see Ness et al., 2013; Bensby et al., 2017; Barbuy et al., 2018). Assuming the metal-poor, low-$\alpha$  stars formed locally, then the radial metallicity gradient would have been $-0.15$ dex kpc*-1* or steeper post gas dilution. This is a much steeper gradient than is measured today. If we use the simple equation

[TABLE]

assuming the present-day metallicity gradient and a difference in [M/H] of 0.9 dex, then the metal-poor, low-$\alpha$  stars currently found in the solar neighborhood must have formed at least 15 kpc farther out than the metal-rich stars.

However, if the inner disc gas was not strongly diluted by the merger event, but continued to form stars from gas enriched mainly by in situ stars, then the metallicity gradient would have been enhanced by the gas accretion and the metal-poor, low-$\alpha$  stars could have migrated a shorter distance. Minchev et al. (2018) find that the radial gradient was likely $-0.15$ dex kpc*-1* at the earliest times and flattened with time, suggesting a gradient of approximately $-0.1$ dex kpc*-1* around the time in question. Such a gradient implies that the metal-poor, low-$\alpha$  stars formed 10 kpc outwards of the metal-rich stars.

Spitoni et al. (2019) recently presented a chemical evolution model that suggests a two infall model is sufficient to reproduce the [$\alpha$/Fe]-[Fe/H]-age distribution of the solar neighborhood using asteroseismic ages. Figure 7 is qualitatively in agreement, but the mean age of the metal-rich stars in this work is 5-6 Gyr compared to 8-10 Gyr as predicted by Spitoni et al. (2019). The significant age difference between the high-$\alpha$  sequence and the metal-rich stars supports the picture of the metal-rich stars having migrated versus having formed locally during the first of two main epochs on star formation. From our data, the latter would require a very extended initial star formation period before the gas infall event. More sophisticated comparisons with GCE models are beyond the scope of this work, but will be addressed in future work.

Figure 7 is in general agreement with similar figures using LAMOST data in Xiang et al. (2017) and Wu et al. (2018). Although these [$\alpha$/M] vs [M/H] distribution are more extended, the age evolution is quite smooth. In Xiang et al. (2017) the high-$\alpha$  sequence transitions from uniformly old into intermediate ages at a lower [M/H] than in Figure 7. Wu et al. (2018) use asteroseismic ages and find the old ages are present until solar [M/H], in agreement with the present work. Ness et al. (2016) and Ho et al. (2017) use [C/N]-based ages produced by the Cannon. Ness et al. (2016) do not find that the high [M/H] stars are older than the solar abundance stars. Ho et al. (2017) do not extend past [M/H] of 0.2, but is consistent with the present results around solar [M/H].

4 Conclusions

By combining the APOGEE spectroscopic survey with the Gaia DR2 astrometric catalogue, the sample of red giant stars for which isochrone matching ages are possible has been vastly increased. The sample presented here contains over 75,000 stars with 5 < $R_{\rm Gal}$  < 13 kpc and 0 < $|z|$  < 2 kpc. The hierarchical modelling method of Feuillet et al. (2016) was used to derive age-abundance trends for [M/H] (MAR) and [$\alpha$/M] (AAR) as a function of spatial location in the disc of the Milky Way. This allows us to examine the spatial variations in age-abundance relations on a disc-wide scale for the first time using ages that do not rely on chemical abundance tracers.

There is significant variation in the MAR through the Milky Way disc, Figure 3, an encouraging result for the potential diagnostic power of such observations to constrain large-scale galaxy simulations and models of chemical evolution. These observations suggest that radial migration has a non-negligible effect in the disc plane. From the metallicity of the youngest stars at each $R_{\rm Gal}$  zone in the plane of the disc, the present-day metallicity gradient is measured to be $-0.059\pm 0.010$ dex kpc*-1*, in agreement with measurements from Cepheids and young field subgiants. The radial metallicity gradient is a key constraint to model of disc evolution (see Stanghellini et al., 2019).

The outer $R_{\rm Gal}$  zones show evidence in the MAR such that young stars are dominant at larger $|z|$  in the outer disc than in the inner disc. This is shown explicitly in Figure 5. The flaring of the outer disc is predicted by models of Milky Way evolution (see Rahimi et al., 2014; Minchev et al., 2015) and observed in large surveys (e.g. Ness et al., 2016; Xiang et al., 2017; Mackereth et al., 2017).

The AAR, Figure 6, also shows evidence of flaring, but is otherwise fairly consistent in all spatial zones. The relative lack of spatial variation in the AAR is also a strong constraint to models of Milky Way disc evolution.

The age trends seen in both [M/H] and [$\alpha$/M] of the solar neighborhood are nicely recovered in Figure 7, which shows the mean age of mono-[M/H], mono-[$\alpha$/M] bins. The high-$\alpha$  sequence is uniformly old with little age evolution until just above solar metallicity. The low-$\alpha$  sequence is younger than the high-$\alpha$  sequence at all [M/H] and shows more significant age evolution. This chemo-age distribution is consistent with previous derivations and models of the Milky Way disc star formation history that find a peak in star formation at around 10 Gyr (forming the high-$\alpha$  sequence) followed by a decline or lull in star formation around 8 Gyr and an extended period of moderate star formation (forming the low-$\alpha$  sequence, e.g. Snaith et al., 2014; Kubryk et al., 2015; Rybizki et al., 2017). The high-[M/H] stars forming the MAR turn over feature are more likely to be an extension of the high-$\alpha$  sequence, perhaps resulting from a continuation of early star formation in the inner disc, rather than an evolution of the low-$\alpha$  sequence.

We plan to use these observations, as well as age-abundance relations of individual elements, to constrain models of chemical evolution and Galactic simulations. Extended disc coverage is crucial for more meaningful comparisons to galaxy scale models.

Acknowledgements

We thank the anonymous referee for useful suggestions. Many thanks to Ted Mackereth and Ása Skúladóttir for helpful discussions.

DKF and KL acknowledge funds from the Alexander von Humboldt Foundation in the framework of the Sofia Kovalevskaja Award endowed by the Federal Ministry of Education and Research. KL also acknowledges funds from the Swedish Research Council (Grant nr. 2015-00415_3) and Marie Sklodowska Curie Actions (Cofund Project INCA 600398). DAGH acknowledges support from the State Research Agency (AEI) of the Ministry of Science, Innovation and Universities (MCIU) and the European Regional Development Fund (FEDER) under grant AYA2017-88254-P. ARL acknowledges financial support provided in Chile by Comisión Nacional de Investigación Científica y Tecnológica (CONICYT) through the FONDECYT project 1170476 and by the QUIMAL project 130001

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-IV acknowledges support and resources from the Center for High-Performance Computing at the University of Utah. The SDSS web site is www.sdss.org.

SDSS-IV 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, the Chilean Participation Group, the French Participation Group, Harvard-Smithsonian Center for Astrophysics, Instituto de Astrofísica de Canarias, The Johns Hopkins University, Kavli Institute for the Physics and Mathematics of the Universe (IPMU) / University of Tokyo, 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. Collaboration Overview Start Guide Affiliate Institutions Key People in SDSS Collaboration Council Committee on Inclusiveness Architects Survey Science Teams and Working Groups Publication Policy How to Cite SDSS External Collaborator Policy