The age-metallicity structure of the Milky Way disk

TL;DR

This study dissects the Milky Way disk into mono-age and mono-[Fe/H] populations using SDSS-APOGEE data, revealing how different stellar populations' structures and contributions vary with age and chemical composition, informing galaxy evolution models.

Contribution

It provides a detailed age- and chemistry-based structural analysis of the Milky Way disk, highlighting the effects of disk heating, radial migration, and population thickness differences.

Findings

Low [$oldsymbol{ m f alpha}$/Fe] populations fit broken exponential profiles.

High [$oldsymbol{ m f alpha}$/Fe] populations are well described by single exponentials.

High [$oldsymbol{ m f alpha}$/Fe] populations are thicker and dominate old stellar mass.

Abstract

The measurement of the structure of stellar populations in the Milky Way disk places fundamental constraints on models of galaxy formation and evolution. Previously, the disk's structure has been studied in terms of populations defined geometrically and/or chemically, but a decomposition based on stellar ages provides a more direct connection to the history of the disk, and stronger constraint on theory. Here, we use positions, abundances and ages for 31,244 red giant branch stars from the SDSS-APOGEE survey, spanning $3 < R_{\mathrm{gc}} < 15$ kpc, to dissect the disk into mono-age and mono-[Fe/H] populations at low and high [$\alpha$/Fe]. For each population, with $\Delta \mathrm{age} < 2$ Gyr and $\Delta \mathrm{[Fe/H]} < 0.1$ dex, we measure the structure and surface-mass density contribution. We find that low [$\alpha$/Fe] mono-age populations are fit well by a broken exponential, which increases to a peak radius and decreases thereafter. We show that this profile becomes broader with age, interpreted here as a new signal of disk heating and radial migration. High [$\alpha$/Fe] populations are well fit as single exponentials within the radial range considered, with an average scale length of $1.9\pm 0.1$ kpc. We find that the relative contribution of high to low [$\alpha$/Fe] populations at $R_0$ is $f_\Sigma = 18\% \pm 5\%$; high [$\alpha$/Fe] contributes most of the mass at old ages, and low [$\alpha$/Fe] at young ages. The low and high [$\alpha$/Fe] populations overlap in age at intermediate [Fe/H], although both contribute mass at $R_{0}$ across the full range of [Fe/H]. The mass weighted scale height $h_Z$ distribution is a smoothly declining exponential function. High [$\alpha$/Fe] populations are thicker than low [$\alpha$/Fe], and the average $h_Z$ increases steadily with age, between 200 and 600 pc.

Figures

Equations

$$\begin{gathered} \lambda (O|\theta )={\nu }_{\ast }(X,Y,Z|\theta )\times |J(X,Y,Z;l,b,D)| \\ \times \rho (H,[J-K_S{]}_0,[Fe/H]|X,Y,Z)\times S(l,b,H) \end{gathered}$$

$$\ln \mathcal{L}(\theta )=\sum _i\left[\ln {\nu }_{\ast }(X_i,Y_i,Z_i|\theta )-\ln \int dO\lambda (O|\theta )\right]$$

$$\begin{gathered} \int dO\lambda (O|\theta )=\sum _{\text{fields}}{\Omega }_f\int dD{D}^2{\nu }_{\ast }([X,Y,Z](D,\text{field})|\theta ) \\ \times \mathfrak{S}(\text{field},D) \end{gathered}$$

$$\begin{gathered} \mathfrak{S}(\text{field},D)=\sum _{k}S(\text{field},k) \\ \int d{M}_H\frac{{\Omega }_k(H_{\text{[min,max]},k},M_H,A_H[l,b,D],D)}{{\Omega }_f}. \end{gathered}$$

$$\begin{gathered} {\Omega }_k(H_{\text{[min,max]},k},M_H,A_H[l,b,D],D)= \\ \Omega (H_{\text{min},k}-[M_H-\mu (D)]<A_H(l,b,D)<H_{\text{max},k}-[M_H-\mu (D)]) \end{gathered}$$

$${\nu }_{\ast }(R,\varphi ,Z)=\Sigma (R)\zeta (Z|R)\text{where }\int dZ\zeta (Z|R)=1.$$

$$\ln \Sigma (R)\propto \{\begin{array}{ll}-h_{R,\text{in}}^{-1}(R-R_0)& \text{where }R\le R_{\text{peak}}\\ -h_{R,\text{out}}^{-1}(R-R_0)& \text{where }R>R_{\text{peak}}\end{array}$$

$$\ln \zeta (Z|R)\propto h_Z^{-1}\exp (R_{\mathrm{flare}}^{-1}[R-R_0])|Z|-\ln h_Z(R).$$

$$N_{R_0}=\frac{N_{\ast ,\text{observed}}}{\int dO\lambda (O|\theta )}$$

$${\Sigma }_{R_0}(\mathrm{age},[Fe/H])=N_{R_0}\frac{⟨M_{\text{RGB}}⟩(\mathrm{age},[Fe/H])}{\omega (\mathrm{age},[Fe/H])}$$

Full text

††pubyear: 2017††pagerange: The age-metallicity structure of the Milky Way disk –B

1 Introduction

The understanding of the present day spatial, kinematic and chemical configuration of the stars of the Milky Way disk is a cornerstone of Galactic archaeology, placing key constraints on models of galaxy disk formation and evolution. Much of our understanding of the time evolution of galaxy disks like that of the Milky Way has arisen from studies which match galaxies of a given stellar mass at $z=0$ to their progenitors at higher $z$ (and therefore, lookback time, e.g. van Dokkum et al., 2013; Papovich et al., 2015; Huertas-Company et al., 2016). However, the Sun’s position in the Milky Way presents a high-fidelity insight into the structure of a Galactic disk on a star by star basis, which has provided a great many insights into the problem (e.g. Eggen et al., 1962; Edvardsson et al., 1993; Haywood et al., 2013). Data for large numbers of disk stars over a wide range of Galactocentric distances, including positions, chemical abundances and stellar ages are now becoming readily available, due to the advent of modern spectroscopic surveys such as APOGEE (Majewski et al., 2015), Gaia-ESO (Gilmore et al., 2012) and GALAH (Martell et al., 2016), with future instruments aiming to bolster the ESA-Gaia data releases (Gaia Collaboration et al., 2016), such as WEAVE (Dalton et al., 2014) and MOONS (Cirasuolo et al., 2012).

Galaxy disks are commonly considered to have stellar density distributions described by exponential laws of some form (e.g. de Vaucouleurs, 1959; Freeman, 1970; Gilmore & Reid, 1983; Pohlen & Trujillo, 2006), assumed classically as the result of gas collapse with angular momentum conservation (e.g. Fall & Efstathiou, 1980), and more recently, the redistribution or ’scrambling’ of the angular momentum of individual stars (e.g. Elmegreen & Struck, 2013, 2016; Herpich et al., 2016). External disks have relatively well constrained photometric scale lengths (e.g. Fathi et al., 2010), but estimates of that of the Milky Way vary greatly (see Bland-Hawthorn & Gerhard, 2016, for a review), and seem to be in tension with those in external galaxies which are assumed to be similar to the Milky Way. This suggests a discord between internal and external scale length measurement methods, or that the Milky Way has a structure distinct from the bulk of disk galaxies.

In the Milky Way, the definition of the measured population appears to have a great effect on this estimate, with thicker, geometrically defined populations having generally flatter radial profiles (e.g. Jurić et al., 2008) than, for example, the populations enhanced in $\alpha$-element abundances (e.g. Cheng et al., 2012b; Bovy et al., 2012b; Bovy et al., 2016b). Theoretical results have suggested that these geometric thick disks are formed from embedded flaring of co-eval populations (Minchev et al., 2015). The discrepant scale lengths between geometric and abundance selected thick disks was framed most recently by Martig et al. (2016a) as the result of a hitherto unaccounted-for radial age gradient in the disk.

When considered in toto, nearby galaxy disks, and the Milky Way disk, are observed to have a two-component vertical spatial structure, commonly referred to as the geometric ‘thick’ and ‘thin’ disks (Tsikoudi, 1979; Burstein, 1979; Yoshii, 1982; Gilmore & Reid, 1983; Jurić et al., 2008). Classically, the ‘thick’ components have been characterised by older, kinematically hotter stellar populations, enriched in $\mathrm{[}\alpha\mathrm{/Fe]}$, whereas the ‘thin’ populations assume lower, near solar $\mathrm{[}\alpha\mathrm{/Fe]}$ and are kinematically cooler (e.g. Bensby et al., 2005). It has, however, also been posited that these populations are in fact composed of multiple sub-populations that smoothly span this range of properties (e.g. Norris, 1987; Nemec & Nemec, 1991; Bovy et al., 2012c, b; Bovy et al., 2016b), thus, in terms of structural parameters, the disk cannot be characterised as the superposition of two distinct structures with different scale heights (Bovy et al., 2012a). It has been shown that the total vertical stellar spatial distribution resulting from the overlap of such sub-populations is consistent with a double exponential (see, e.g., Figure 14 of Rix & Bovy, 2013). It is difficult to explain the presence of a continuity in structure alongside the discontinuity in chemistry seen in the Milky Way.

The Milky Way’s disk has a complex chemical structure, with a bimodality in $\mathrm{[}\alpha\mathrm{/Fe]}$ seen at fixed $\mathrm{[Fe/H]}$ across many of the observable regions of the disk (Bensby et al., 2003, 2005; Nidever et al., 2014; Hayden et al., 2015). This characteristic is difficult to explain using one-zone Galactic chemical evolution (GCE) models (most recently shown by Andrews et al., 2017), giving rise to attempts to explain it by means other than pure chemical evolution. Examples of such models include the heating of the old disk by high redshift mergers (e.g. Brook et al., 2004; Villalobos & Helmi, 2008; Kazantzidis et al., 2009; Minchev et al., 2013) and the formation of a dual disk by gradual accretion of stars into disk orbits (e.g. Abadi et al., 2003). More recent work has also framed this bimodality as a consequence of discontinuous radial migration of stars in the disk (Toyouchi & Chiba, 2016).

However, such chemical structure can be replicated in part by invoking various Galactic chemical evolution models that do not rely on a ‘one-zone’ approximation (e.g. Chiappini et al., 1997; Portinari & Chiosi, 2000; Andrews et al., 2017; Weinberg et al., 2017). For example, Andrews et al. (2017) showed that a combination of GCE models with varying outflow mass loading parameters and inflow timescales (intended to represent enrichment histories at varying Galactocentric radii) could make a roughly bimodal $\mathrm{[}\alpha\mathrm{/Fe]}$ distribution. The same models were shown to present a good explanation of the APOGEE $\mathrm{[}\alpha\mathrm{/Fe]}$-$\mathrm{[Fe/H]}$ plane by Nidever et al. (2014). A deeper understanding of the connection between spatial structure in $\mathrm{[}\alpha\mathrm{/Fe]}$ and stellar age selected populations in the Milky Way is necessary to link these results.

In this paper, we present the first dissection of radially extended samples of Milky Way disk stars in age, $\mathrm{[Fe/H]}$ and $\mathrm{[}\alpha\mathrm{/Fe]}$. A strong correlation is observed between stellar age and $\mathrm{[}\alpha\mathrm{/Fe]}$ in the solar vicinity (Haywood et al., 2013). On the other hand, but also in the solar vicinity, no correlation is found between age and [Fe/H] (e.g. Edvardsson et al., 1993; Nordström et al., 2004), which may be explained by the occurrence of radial migrations. Thickening of the Galactic disk has been invoked as a consequence of outward stellar radial migration (e.g. Schönrich & Binney, 2009b). However, Minchev et al. (2012) argue that the effect is small and that such migration in fact only makes disks flare by a small amount. Similarly, Bovy et al. (2016b) measured the flaring profile of low $\mathrm{[}\alpha\mathrm{/Fe]}$ stars to only slowly exponentially increase with Galactocentric radius, and suggest that radial migration is likely not a viable mechanism for forming thickened disk components. Flaring has also been shown to arise as a result of satellite infall, which can be a stronger flaring agent than migrations (e.g. Bournaud et al., 2009). The understanding of flaring and its connection to the evolution of the Galactic disk is essential, but as yet incomplete. In this paper, we present new constraints on models of radial migration in the disk by studying its effects on the Milky Way’s mono-age stellar populations.

Many theoretical studies have attempted to understand the observed structure of mono-abundance populations (MAPs) through the use of hydrodynamics and N-body simulations. Few reproduce the observed bimodality in $\mathrm{[}\alpha\mathrm{/Fe]}$ at fixed $\mathrm{[Fe/H]}$, and so an understanding of this has so far proved difficult. However, certain characteristics of the Milky Way $\mathrm{[}\alpha\mathrm{/Fe]}$ distribution are beginning to emerge in the most recent cosmological simulations (Ma et al., 2016). Structurally, the mono-age populations of simulated galaxies show good agreement with the Milky Way (e.g. Stinson et al., 2013; Bird et al., 2013; Martig et al., 2014a, b). More recent work has brought into question the applicability of MAPs as a proxy for mono-age populations (Minchev et al., 2017), showing that, particularly at low $\mathrm{[}\alpha\mathrm{/Fe]}$, MAPs may have significant age spreads due to the differential nature of star formation in the disk. We show in this paper that the structures of mono-age and mono-abundance populations in the Milky Way disk differ, but present complementary insights into the temporal and chemical evolution processes in the disk.

Previous work has studied MAPs in the Milky Way by analysing samples of SEGUE G-dwarfs (Bovy et al., 2012c, b; Bovy et al., 2012a) and APOGEE red-clump (RC) giants (Bovy et al., 2016b). In this work we map the spatial distribution of mono-age, mono-$\mathrm{[Fe/H]}$ populations at low and high $\mathrm{[}\alpha\mathrm{/Fe]}$ using a catalogue containing $\mathrm{[Fe/H]}$ and $\mathrm{[}\alpha\mathrm{/Fe]}$ from the APOGEE survey (Majewski et al., 2015) and ages from Martig et al. (2016b) for 31,244 red giant stars. We complement earlier work by adapting the method developed by Bovy et al. (2016b) for RC stars, to enable its application to the full red giant branch (RGB) sample from APOGEE. Stars in the RGB are better tracers of the underlying stellar population than their RC counterparts because of reduced uncertainties in the stellar evolution models, and because they are generally brighter, and so are observed at greater distances. On the other hand, this means that the method developed by Bovy et al. (2016b) must be adapted to account for the spread in absolute magnitude over the RGB (whereas RC stars can be considered as a near-standard candle). On the basis of these measurements, we establish the local mass-weighted age-$\mathrm{[Fe/H]}$ distribution, showing the contributions from both low and high $\mathrm{[}\alpha\mathrm{/Fe]}$ stellar populations.

In Section 2, we discuss the APOGEE data and the distance and age catalogues used for this work. Section 3 describes the stellar density fitting method, drawing greatly on work by Bovy et al. (2016a); Bovy et al. (2016b). Specifically, we describe the generalities of the maximum likelihood fitting procedure and the calculation of the effective survey selection function for RGB stars in Section 3.1, the adopted parametric stellar density model in Section 3.2, and the method for calculating stellar surface-mass densities in Section 3.3. We present the fits in Section 4, along with the calculated surface-mass density contributions for each mono-age, mono-$\mathrm{[Fe/H]}$ population. In Section 5 we compare our findings to those in the literature and discuss possible scenarios for the formation of the Milky Way’s disk in light of our findings. Section 6 summarises our results and conclusions.

2 Data

2.1 The APOGEE Catalogue

We use data from the twelfth data release (DR12, Alam et al., 2015) of the SDSS-III APOGEE survey (Majewski et al., 2015), a high signal-to-noise ratio (SNR$>100\ \mathrm{pixel}^{-1}$), high resolution ($R\sim 22,500$), spectroscopic survey of over 150,000 Milky Way stars in the near-infrared $H$ Band ($1.5-1.7\mathrm{\mu m}$). Stars were observed during bright time with the APOGEE spectrograph (Wilson et al., 2010) on the 2.5m Sloan Foundation Telescope (Gunn et al., 2006) at Apache Point Observatory. Targets were selected in general from the 2MASS point-source catalog, employing a dereddened $(J-K_{S})_{0}\geq 0.5$ colour cut (in the fields which are of interest here) in up to three apparent $H$ magnitude bins (for a full description of the APOGEE target selection, see Zasowski et al., 2013). Reddening corrections were determined for the colour cut via the Rayleigh-Jeans Colour Excess method (RJCE, Majewski et al., 2011). Corrections are found by applying the method to 2MASS (Skrutskie et al., 2006) and mid-IR data from Spitzer-IRAC GLIMPSE-I, -II, and -3D (Churchwell et al., 2009) when available and from WISE (Wright et al., 2010) otherwise. For this work, we use distance moduli (which make use of the aforementioned reddening corrections) from the Hayden et al. (2015) distance catalogue for DR12 (see Section 2.2).

All APOGEE data products employed in this paper are those output by the standard data reduction and analysis pipeline used for DR12. The data were processed (Nidever et al., 2015), then fed into the APOGEE Stellar Parameters and Chemical Abundances Pipeline (ASPCAP, García Pérez et al., 2016), which makes use of a specifically computed spectral library (Zamora et al., 2015), calculated using a customised $H$-band line-list (Shetrone et al., 2015). Outputs from ASPCAP are analysed, calibrated and tabulated (Holtzman et al., 2015). The output $\mathrm{[}\alpha\mathrm{/Fe]}$ and $\mathrm{[Fe/H]}$ abundances for DR12 have been shown to have a high degree of precision (at least between $4500\lesssim\mathrm{T_{\mathrm{eff}}}\lesssim 5200$ K), such that $\sigma_{\mathrm{[Fe/H]}}=0.05\ \mathrm{dex}$ and $\sigma_{\mathrm{[\alpha/Fe]}}=0.02\ \mathrm{dex}$ (Bovy et al., 2016b). We apply here the same external calibrations to $\mathrm{[}\alpha\mathrm{/Fe]}$ and $\mathrm{[Fe/H]}$ as Bovy et al. (2016b), constant offsets of -0.05 dex and -0.1 dex, respectively. We use the tabulated $\mathrm{[Fe/H]}$ value rather than the globally fit $\mathrm{[M/H]}$ which is included in the table.

We select stars from the DR12 catalogue which were targeted as part of the main disk survey (i.e were subject to the $(J-K_{S})_{0}\geq 0.5$ cut), have reliably measured abundances (i.e. no warning or error bits set in the ASPCAPFLAG field) and have a well defined distance modulus and age measurement (see Sections 2.2 and 2.3). We apply a secondary cut at $1.8<\log{g}<3.0$ to restrict the sample to stars on the red giant branch (RGB), removing most contaminating dwarfs, and very evolved stars near the tip of the RGB. The high end of our $\log{g}$ cut is more conservative than other studies in this regime, however we find this gives the best agreement between the data and stellar evolution models, without significantly reducing the sample size or introducing unwanted bias. These cuts give a final sample of 31,244 stars, spanning $4200\lesssim\mathrm{T_{\mathrm{eff}}}\lesssim 5050\ \mathrm{K}$ for which the effective survey selection function can be reconstructed.

We further divide the sample into low and high $\mathrm{[}\alpha\mathrm{/Fe]}$ sub-samples, as it has been shown in previous work that the two populations have quite different structural parameters (Bovy et al., 2016b), and as such it makes sense to fit their mono-age sub-populations separately. We separate visually the low and high $\mathrm{[}\alpha\mathrm{/Fe]}$ populations, leaving a gap between the two samples of 0.05 dex in $\mathrm{[}\alpha\mathrm{/Fe]}$ at each $\mathrm{[Fe/H]}$ (our separation is shown in Figure 1), to minimise contamination between the subsamples, particularly at the high $\mathrm{[Fe/H]}$ end, where the two populations partially overlap. The final density fits are performed on finer bins in age and $\mathrm{[Fe/H]}$, which we define in Section 2.3. As the adopted separation in $\mathrm{[}\alpha\mathrm{/Fe]}$ removes 6532 stars from the full count, when calculating the surface-mass density contributions from the stellar number counts in each age-$\mathrm{[Fe/H]}$ bin, we remove the separation in $\mathrm{[}\alpha\mathrm{/Fe]}$, using the star counts as if the populations were separated along the midpoint of the division (as shown by the dot-dashed line in Figure 1).

Our method, discussed in Section 3, corrects for selection effects induced by interstellar extinction (in addition to the RJCE reddening corrections) using 3D dust maps for the Milky Way derived by Marshall et al. (2006) for the inner disk plane, combined with those for a large majority of the APOGEE footprint by (Green et al., 2015), adopting conversions $A_{H}/A_{K_{S}}=1.48$ and $A_{H}/E(B-V)=0.46$ (Schlafly & Finkbeiner, 2011; Yuan et al., 2013). Fields with no dust data (of which there are $\sim 10$) are removed from the analysis. Bovy et al. (2016b) discuss the relative merits and limitations of these dust maps as opposed to others which are available, and determine that this combination of dust maps provides the best density fits.

2.2 Distance estimates

We use distance estimates from Hayden et al. (2014) (But see also Hayden et al., 2015, for further description). Distances are estimated by computing the probability distribution function (PDF) of all distance moduli to a given star using a Bayesian method applied to the spectroscopic and photometric parameters from the DR12 catalogue and the PARSEC isochrones (Bressan et al., 2012). The distance estimates are found to have accuracy at the $15-20\%$ level, upon comparison with cluster members of well known distance observed by APOGEE.

We use the median of the posterior PDF (which is given in the output catalogue) as the estimate for the distance modulus, and compute the Galactocentric cylindrical coordinates, $R,\phi\ \rm{and}$ $Z$ for each star using the $l,b$ coordinates provided in the APOGEE-DR12 catalogue. The spatial distribution of the two $\mathrm{[}\alpha\mathrm{/Fe]}$ sub-samples is shown in Figure 2. We perform a simple cross match between our sample and the APOGEE red clump (RC) value added catalogue (VAC) (Bovy et al., 2014) and plot the red-clump derived distance ($D_{RC}$) against the estimate from Hayden et al. (2014) ($D_{MH}$) in Figure 3. The red clump catalogue has precise distance estimates, which can be determined due to the red-clump having a near constant absolute magnitude. The majority of the Hayden et al. (2014) distances compare well to the RC distances, but there are notable differences. The Hayden et al. (2014) distances can be underestimated by as much as 50%, and we find that $\sim 20\%$ of our sample have distances underestimated by more than 10%. Our density fitting method is insensitive to uncertainties on the distances at these scales, as the scale of any variations in the density distribution can be assumed to be far greater than the distance uncertainties.

Figure 3 also shows that there is a systematic offset between the RC and Hayden et al. (2014) distances of the order $\sim 5\%$ across the full range of distances. We find that adopting this offset as a correction to the distances makes little impact on the final results, merely broadening fitted density profiles slightly, spreading the star counts over a wider Galactocentric distance, meaning that the final stellar surface-mass density estimates are unchanged.

2.3 Age estimates

We use age estimates for APOGEE DR12 catalogued by Martig et al. (2016b), who derive an empirical model for the $\mathrm{[C/N]}-M_{*}$ relation using asteroseismic masses from Kepler and abundances from APOGEE for their overlapping samples (APOKASC Pinsonneault et al., 2014). Masses are predicted for DR12 stars which meet quality and stellar parameter criteria outlined in Martig et al. (2016b), and ages are estimated from that mass using the PARSEC isochrones with the nearest metallicity to that of the given star. Martig et al. (2016b) use this empirical relation to build a model which predicts mass and age as a function of $[\text{[M/H], [C/M], [N/M], [(C+N)/M],}\log{g},T_{\mathrm{eff}}]$. It is important to note here that Martig et al. (2016b) derive a model and fit for the ages in DR12 using the uncalibrated, raw stellar parameters, found in the FPARAM arrays in the APOGEE catalogue. This is difficult to account for when using the age catalogue alongside the calibrated parameters, and must be borne in mind in future comparisons of this work with models and observational results. In addition to this, Martig et al. (2016b) also mention that care should be taken when applying these ages to regions of the Milky Way where the chemical evolution may have been complex (e.g. the Bulge/Bar region). However, in their Figure 12, they compare the [C/N] ratio as a function of [M/H] in a sample of pre-dredge-up giants in the inner and outer disk, showing that the shapes of the distributions are similar. This suggests that differences in chemical evolution do not affect the $\mathrm{[C/N]}$-age relation within a wide range of galactocentric distances. Therefore, the assumption that it is safe to adopt the Martig et al. (2016b) ages over the extent of the disk covered by our sample is robust, regardless of the fact that they are trained on the Kepler sample, which is limited in its spatial extent.

Although individual uncertainties on ages are not given in the catalogue, Martig et al. (2016b) state that the model predicts ages with r.m.s errors of $\sim 40\%$. Although uncertainties are potentially very large at high age, our sample is binned with $\Delta\text{age}=2$ Gyr in order to gauge general trends with age. It should be understood that such trends are smoothed by the age uncertainties, particularly at high age, and detailed comparisons to models should take the age uncertainty into account. We discuss the effect of these uncertainties on our recovered trends with age in Appendix B, showing that our methodology can still reliably recover such trends, even though mixing between bins may be present in the data.

Figure 11 of Martig et al. (2016b) shows that there is a significant bias in the ages returned by the model, such that ages are underpredicted at high age when compared to the training set. For this reason, we fit for and apply a correction to the catalogued ages before performing the density fitting. Using Table 1 of Martig et al. (2016b), we perform a non-parametric lowess fit to the predicted age–true age distribution. This fit is then used to derive age corrections as a function of predicted age. We show the fitted correction in Figure 4 in both predicted vs. true age space and also in $\Delta$ age against the predicted age. The correction as a function of predicted age is then applied to each of the ages in the DR12 catalogue. In all further analysis, we refer only to the corrected ages. The main effect of this correction is to make the high-$\mathrm{[}\alpha\mathrm{/Fe]}$ stars older. Consequently, our surface-mass density estimates (presented in Section 4.3) become more conservative, as the mass contribution per star in older bins is lower (as discussed in Section 5.1). Martig et al. (2016b) comment on the bias in the context that the ages returned for high $\mathrm{[}\alpha\mathrm{/Fe]}$ stars appear younger than previous estimates (from Haywood et al., 2013; Bensby et al., 2014; Bergemann et al., 2014). Our correction brings these data more in line with those estimates.

It is also important to account for all other cuts made by Martig et al. (2016b) on the stellar parameters in the APOGEE catalogue, outlined in full at the beginning of their Section 6.2. While we account for cuts in $T_{\mathrm{eff}}$ and $\log{g}$ by applying the same cuts to the isochrone grid when calculating surface-mass density contributions, it is not possible to properly account for cuts made on the stellar abundances in this way. We find that 9,041 stars are removed from the 40,285 star catalogue (those with distances, after the $\log{g}$cut mentioned above) by these abundance cuts, to give our final catalogue size of 31,244. This means that $\sim 25\%$ of star counts are missing from the age catalogue, and therefore unaccounted for by our analysis. With no robust method for determining the age distribution of these missing stars, we are forced to make the assumption that these star counts can be added uniformly to each age-$\mathrm{[Fe/H]}$ bin. We make this correction by simply increasing the counts in each bin by $25\%$ when calculating the surface-mass density. This correction simply increases the final surface-mass density values systematically by $25\%$.

Another important consideration when using this set of ages regards stars whose chemical compositions are such that the ages fit from the model (after making the corrections) would be higher than 13 Gyr and left out of our analysis. As the age estimates are computed based on the surface parameters and abundances of the stars using a fitting function, many stars with strongly outlying abundances and parameters can be assigned ages which are greater than 13 Gyr. We find that the 3020 such stars (in the final sample) have an average [C/Fe] which is lower than the general sample, and an average [N/Fe] which is enhanced with respect to the stars with reliably measured ages. We also find that at fixed $\mathrm{[Fe/H]}$ and $\log{g}$ these stars have warmer $T_{\mathrm{eff}}$. While these properties are expected given their age measurements, we believe that there is a distinct possibility of some peculiarity of these stars. For example, if such stars were early AGB stars (having gone through the second dredge-up, reducing the surface C abundance, e.g. Boothroyd & Sackmann, 1999), which had been fit as RGB stars, their actual age may be considerably younger, and their counts missed in the younger bins (where mass contribution per star is higher). As the nature of these stars is debatable and a correction cannot be made confidently before carrying out the full analysis, we regard the missing counts from these stars as a contribution to the systematic error budget, which we discuss fully in Section 5.1.

We demonstrate the adopted binning in (age,[Fe/H]) space in Figure 5, showing also the number of stars which fall in each bin and the general distribution of stars in (age,[Fe/H]) space. There is a notable separation in age between the high and low $\mathrm{[}\alpha\mathrm{/Fe]}$ subsamples.

3 Method

In this section we describe the method for fitting the underlying number density of stars in the Milky Way from APOGEE observations, which we represent here as $\nu_{*}(X,Y,Z|\theta)$, in units of stars kpc*-3*. The calculation of this quantity requires allowances to be made for the survey selection function, which is non-trivial due to the presence of inhomogeneous dust extinction along lines of sight observed by APOGEE, the target selection invoking different $H$ magnitude limits, and the use of RGB stars as a tracer, which cannot be considered as standard candles. The quantity which we are ultimately interested in is the surface-mass density of stars at the solar radius, $\Sigma_{R_{0}}$, in units of $\mathrm{M_{\odot}}\ \mathrm{pc^{-2}}$, which we infer from the number of stars in the APOGEE sample as a function of position. We describe the method for this calculation in Section 3.3. Our methodology consists of an adaptation of that used by Bovy et al. (2016b), employing a modified version of their publicly available code111Available at https://github.com/jobovy/apogee-maps. Although the general method is identical, we describe again the key components for clarity and completeness.

As some readers may find it unnecessary to read in full the details of the methodology (which are described in the following sections), we summarise the procedure as follows:

• •

We fit parametric density models to the APOGEE star counts using a maximum likelihood fitting procedure, based on the assumption that star counts are well modelled as an inhomogeneous Poisson point process. The density models which we assume throughout the paper are described by radially broken exponentials, with scale lengths $h_{R,\mathrm{[in,out]}}$ either side of a break radius $R_{\mathrm{peak}}$ (where $h_{R,\mathrm{in}}$ denotes the scale length of the inner profile and vice versa), and a vertical distribution which is a single exponential with scale height $h_{Z}$, which is modified as a function of $R$ by an exponential flaring term with scale length $R_{\mathrm{flare}}$. We show that, in general, if the density is better fit by a single exponential, it is recovered as so by our procedure.

• •

We obtain a best fit density model for every bin in age and $\mathrm{[Fe/H]}$, at high and low $\mathrm{[}\alpha\mathrm{/Fe]}$. This best fit model is then used to initiate an MCMC sampling of the posterior PDF. We then use the median and standard deviation of one dimensional projections of the MCMC chain as our adopted parameter values and uncertainties.

• •

As the fitting procedure does not fit for the normalisation of the density $N_{R_{0}}$, the number surface density of stars at the solar radius in stars pc*-2*, we calculate this value by comparing the observed number of stars in each bin to that which would be observed in APOGEE for the fitted density model if $N_{R_{0}}=1$ star pc*-2*. We then convert $N_{R_{0}}$ for each bin into the surface-mass density in visible stars at the solar radius $\Sigma_{R_{0}}$ by converting the mass in RGB stars observed to the total mass using stellar evolution models.

Readers can then pick up the results in Section 4.

3.1 Density Fitting Procedure

We first fit for the number density of stars for each sub-population as defined by Figure 5. The following discussion describes the general procedure used for fitting density models with a generic set of parameters $\theta$. The actual stellar number density model adopted is discussed in Section 3.2. Bovy et al. (2016b); Bovy et al. (2012b); Rix & Bovy (2013) have shown that the observed rate of stars as a function of position, magnitude, colour and metallicity can be modelled as an inhomogeneous Poisson point process. Stars are distributed in the space defined by $O=[l,b,D,H,\mathrm{[J-K_{S}]_{0}},\mathrm{[Fe/H]}]$ – position, magnitude, colour and metallicity – with an expected rate $\lambda(O|\theta)$ (which has units of stars per arbitrary volume in O), parameterised by a set of parameters $\theta$ (which are, in this particular case, the parameters describing an adopted density profile). This rate function is written fully as

[TABLE]

where $\nu_{*}(X,Y,Z|\theta)$ is the quantity we aim to estimate, which is defined as the stellar number density in rectangular coordinates, in units of stars kpc*-3*. $|J(X,Y,Z;l,b,D)|$ is the Jacobian of the transformation from rectangular $(X,Y,Z)$ to Galactic $(l,b,D)$ coordinates and $\rho(H,[J-K_{S}]_{0},\mathrm{[Fe/H]}|X,Y,Z)$ denotes the density of stars in magnitude, colour and metallicity space given a spatial position $(X,Y,Z)$, in units of stars per arbitrary volume in magnitude, colour and metallicity space. $S(l,b,H)$ is the survey selection function (the fraction of stars observed in the survey) which includes dust extinction effects, which we discuss in the following. When expressed in this way, fitting the density model parameters $\theta$ becomes a maximum likelihood problem.

The likelihood is a sum over all data-points considered in a given age-$\mathrm{[Fe/H]}$ bin, and gives the likelihood of the parameters $\theta$ given the data. For this application, it is written as

[TABLE]

where second term on the right hand side of the equation, $\int dO\lambda(O|\theta)$, describes the effective volume of the survey. we drop the other factors in the rate in Equation (1) in the argument of the logarithm, because the other factors do not depend on the model parameters $\theta$. The effective volume is independent of the data-point considered, and is an intrinsic property of the survey for a given $\theta$. It provides the normalisation for the rate likelihood, and is non-trivial to evaluate due to the presence of patchy dust extinction along lines of sight in the survey.

The effective volume is written generally as

[TABLE]

which is a sum over all APOGEE fields, where $\Omega_{f}$ is the solid angle of the field considered. The integrand $\nu_{*}([X,Y,Z](D,\text{field})|\theta)$ is the density at each point along a line of sight, assumed to be constant over the angular size of the field. $\mathfrak{S}(\text{field},D)$ represents the effective survey selection function, which is given by the integration of the survey selection function over the area of the field and is written, in this case, as

[TABLE]

This is a sum over the apparent magnitude bins, $k$, in the APOGEE target selection, with the integral representing the fractional area of the APOGEE field where stars are observable, given the distance modulus and extinction at a given position. The term describing this area is $\Omega_{k}$, which is the observable area of the field at a given distance and absolute magnitude, written as

[TABLE]

where $H_{\text{[min,max]},k}$ denotes the minimum and maximum $H$ for an apparent magnitude bin $k$ in the APOGEE target selection and $\mu(D)$ is the distance modulus at $D$. $A_{H}(l,b,D)$ is the $H$ band extinction at a given position, which we obtain from the 3D dust maps described in Section 2.1. This area is integrated (in Equation (4)) over the full absolute $H$-band magnitude, $M_{H}$, distribution in an (age,$\mathrm{[Fe/H]}$) bin. We find the $M_{H}$ distribution for each (age,$\mathrm{[Fe/H]}$) bin using the PARSEC isochrones (Bressan et al., 2012) within that bin, weighted with a Chabrier (2003) IMF. We apply the same cuts in $\log{g}$ and $(J-K_{S})_{0}$ colour to the isochrone points as are imposed on the data, and perform a Monte Carlo integration using the resulting $M_{H}$ distribution to evaluate the integral in Equation (4). S(field,k) in Equation (4) denotes the ’raw’ APOGEE selection function, which gives the fraction of the stars in the photometric catalogue that were observed spectroscopically (see Zasowski et al., 2013, for details). This number is constant within an apparent magnitude bin and within an APOGEE field, which is why $S$ is cast as a function of field and magnitude bin in Equation (4). The values of S(field,k) (and $\mathfrak{S}(\text{field},D)$) are evaluated using the apogee python package 222Available at https://github.com/jobovy/apogee.

We evaluate $\mathfrak{S}(\text{field},D)$ on a grid of distances for each APOGEE field for simple computation of $\int dO\lambda(O|\theta)$. We then optimise the likelihood function in Equation (2) for a given density model and data-set using a downhill-simplex algorithm, to obtain the best fitting set of parameters $\theta$. A Markov Chain Monte Carlo (MCMC) sampling of the posterior PDF is then initiated using this optimal solution. This is implemented with an affine-invariant ensemble MCMC sampler (Goodman & Weare, 2010; Foreman-Mackey et al., 2013). All parameter values and associated uncertainties for individiual (age,$\mathrm{[Fe/H]}$) bins which are reported in the following sections represent the median and standard deviation $\sigma$, respectively, of one dimensional projections of the MCMC chain.

3.2 Adopted stellar number density models

It was shown in Bovy et al. (2016b) that density profiles of MAPs are well represented by axisymmetric profiles that can be written as

[TABLE]

Furthermore, the exact form of the best fitting profile is that of a radially broken exponential, with a vertical profile that is an exponential with a scale height which varies exponentially with R (a flaring profile), such that

[TABLE]

and

[TABLE]

$R_{0}$ denotes the solar radius, which we assume here to be 8 kpc. This number only sets the radius at which the profiles are normalised, and so does not have any effect on the fitting procedure. We use the same general set of density profiles to describe the mono-age, mono-metallicity populations which are studied here. We bin in $\mathrm{[Fe/H]}$ to account for the observed $\mathrm{[Fe/H]}$ spread at fixed age (e.g. Edvardsson et al., 1993, and our Figure 5). Bovy et al. (2016b) also showed that when mock data were fit using the procedure in Section 3.1 and the density profile above, the input parameters were always recovered within acceptable uncertainty ranges. In particular, mock data generated from a single exponential profile was still recovered as such (i.e. with $R_{\mathrm{peak}}=0$) even when fit assuming a broken exponential profile.

We also note here that our sample is not limited to stars which are members of any specific Galactic component, and as such, may include small numbers of halo stars in the very high $\mathrm{[}\alpha\mathrm{/Fe]}$ and low $\mathrm{[Fe/H]}$ regimes. However, our fitting procedure is agnostic to these contaminants, which would only cause the fits to have larger uncertainty (from the MCMC exploration) about the best fit from the dominant population in a given bin.

3.3 Stellar surface-mass densities

We compute the surface-mass density in visible stars for each of our age and $\mathrm{[Fe/H]}$ populations using the method originally outlined in Bovy et al. (2012a). As our fitting procedure does not fit for the normalisation of the density (we normalise to a surface density of 1 at $R_{0}$), we first compute the normalisation $N_{R_{0}}$, which represents the number density of stars at the solar radius in units of stars pc*-2* in an (age,$\mathrm{[Fe/H]}$) bin. $N_{R_{0}}$ is given by the relation

[TABLE]

where $N_{*,\text{observed}}$ is the number of stars observed in the survey for a given (age,$\mathrm{[Fe/H]}$) and $\int dO\lambda(O|\theta)$ is the usual definition of the effective volume (given by Equation (3)) for a given set of parameters $\theta$, found using the method in Section 3.1.

We find the contribution to stellar surface-mass density by first multiplying $N_{R_{0}}$ by the average mass of a red giant star in the same range of age and $\mathrm{[Fe/H]}$, given the selection criteria on $\log{g}$ given in Section 2.1 (which picks out the RGB) and $(J-K_{S})_{0}\geq 0.5$ (given that we only use fields in which this cut was applied). We then correct this value to represent the total stellar population by dividing by the fractional contribution of the red giants to the total underlying population. These values are found using PARSEC isochrones (Bressan et al., 2012), weighted with a Log-normal Chabrier (2001) IMF, as described in the calculation of the effective volume in Section 3.1. This then leads us to a stellar surface-mass density $\Sigma_{R_{0}}$ as a function of age and $\mathrm{[Fe/H]}$. This conversion can be expressed as

[TABLE]

Where $\langle M_{\text{RGB}}\rangle(\mathrm{age,[Fe/H]})$ is the mean stellar mass in an (age,$\mathrm{[Fe/H]}$) bin, and $\omega(\mathrm{age,[Fe/H]})$ is the fraction of stars in the total stellar population in an (age,$\mathrm{[Fe/H]}$) bin which are within the $\log{g}$ and $(J-K_{S})_{0}$ cuts in APOGEE. The stellar surface-mass density contributions of each bin can then be summed to give the total stellar surface-mass density at the solar radius $\Sigma_{R_{0},\text{tot}}$.

Our final surface-mass density estimate is strongly dependent on the conversion factors in the above equations, the average RGB star mass, $\langle M_{\mathrm{RGB}}\rangle$, and the fractional contribution from giants, $\omega$. We find that the average giant masses in our range of ages and $\mathrm{[Fe/H]}$ span $0.9\lesssim\langle M_{\mathrm{RGB}}\rangle\lesssim 2.1\mathrm{M_{\odot}}$. The most metal poor and oldest populations have the lowest average mass, and the youngest, most metal rich populations have the highest. The fractional contribution from giants in this regime ranges between $0.002\lesssim\omega\lesssim 0.02$. The oldest and most metal poor populations have the least giants, whereas the youngest, metal rich populations have the most. These values appear to sit well with recent inventories of the solar neighbourhood, which suggest giants should make up of the order of a few percent of the mass (McKee et al., 2015). We discuss the potential systematics introduced by the use of stellar evolution models in Section 5.1.

4 Results

We now present results from the density fitting procedure, and the subsequent calculation of the surface-mass density contribution of each mono-age, mono-$\mathrm{[Fe/H]}$ population in Figure 5. Density fitting is performed on all populations, but we only display the fits for populations with $>30$ stars, as as data below this level become too noisy to render reliable fits. Although the remaining fits can be noisy when star counts are near this limit, this is reflected in the error analysis arising from the MCMC exploration of the posterior PDF of the fitted parameters. We refer the reader to Appendix A for a comparison between the data and the fitted models for each mono-age, mono-$\mathrm{[Fe/H]}$ bin, and a qualitative discussion regarding the rationale behind the decision to discuss fits to only the broken exponential density profile.

4.1 The radial profile of mono-age, mono-$\mathrm{[Fe/H]}$ populations

We first show the fits to the surface density in the low and high $\mathrm{[}\alpha\mathrm{/Fe]}$ sub-samples in Figure 6. We display fits for all age and $\mathrm{[Fe/H]}$ bins with $>30$ stars. By shading the profiles by their surface-mass density contribution (as shown in Section 4.3), we intend to draw the eye to the profiles which contribute most to the mass of the Milky Way disk. We defer a discussion of the individual mass contributions of each bin to Section 4.3, concentrating in this section on trends in the shapes of the density profiles.

Although fit with the broken exponential, high $\mathrm{[}\alpha\mathrm{/Fe]}$ profiles are generally better described by near-single exponentials, either showing no break in the radial range, or being fit by a profile with a break at low significance (i.e. a single line could be drawn through the coloured band). Many of the outer profiles (after $R_{\mathrm{peak}}$) in the high $\mathrm{[}\alpha\mathrm{/Fe]}$ sub-sample appear to have a similar slope, suggesting that they may all be represented by the same exponential. The mean outer scale length for the high $\mathrm{[}\alpha\mathrm{/Fe]}$ populations is $h_{R,\text{out}}=1.9\pm 0.1$ kpc. The picture is noticeably different in the low $\mathrm{[}\alpha\mathrm{/Fe]}$ sub-sample, with profiles showing clear breaks, at well defined radii. Any trends in break radius in this regime are determined with high significance. Low $\mathrm{[}\alpha\mathrm{/Fe]}$ profiles have a density which increases with radius out to the break radius, and declines outward of this radius. We do not constrain the fits to behave in this way, and this indicates that mono-age, mono-$\mathrm{[Fe/H]}$ populations at low $\mathrm{[Fe/H]}$ are shaped approximately as donut-like annuli. The variation of the break radius then represents the moving peak of stellar density as a function of age and $\mathrm{[Fe/H]}$. Concentrating on the bins youngest bin ($1<\mathrm{age}<3$ Gyr), the break radius is a declining function of metallicity, moving between $R_{\mathrm{peak}}=10$ kpc at $-0.6<\mathrm{[Fe/H]}<-0.5\ \mathrm{dex}$ down to $R_{\mathrm{peak}}<8$ kpc at $0.1<\mathrm{[Fe/H]}<0.2\ \mathrm{dex}$. This trend is also present in older bins but with decreased amplitude. In a fixed $\mathrm{[Fe/H]}$ bin, $R_{\mathrm{peak}}$ appears to remain roughly constant (within $\sim 1$ kpc) at ages between 1 and 6 Gyr. At ages older than this $R_{\mathrm{peak}}$ varies in unexpected ways, but there is much less mass contribution from these populations, and we attribute much of this behaviour to noise from the narrow age bins.

On the other hand, the low $\mathrm{[}\alpha\mathrm{/Fe]}$ profiles change shape (either side of $R_{\mathrm{peak}}$) with age in a fixed $\mathrm{[Fe/H]}$ bin. The youngest populations show a sharp peak, with a steep increase and decline either side of $R_{\mathrm{peak}}$. As populations grow older, the profile broadens significantly, becoming almost flat in the lowest $\mathrm{[Fe/H]}$ bins. We show this behaviour by finding the inverse of the difference between the inverse outer and inner scale length333Taking a ratio of the sum of the density at fixed $\Delta R$ either side of $R_{\mathrm{peak}}$ to that at $R_{\mathrm{peak}}$ would give some measure of width. Then, assuming $\Delta R<<h_{R,[\text{in,out}]}$, a Taylor expansion of this ratio $\sim\Delta R(h_{R,\text{out}}^{-1}-h_{R,\text{in}}^{-1})$. We then plot the inverse of this factor such that it increases for broader profiles., such that a low value denotes a sharper peak, whereas a broader profile has a higher value. We show how this value changes with age for the low $\mathrm{[}\alpha\mathrm{/Fe]}$ populations in Figure 8. The peak is sharpest in the younger populations, and becomes broader with age. Old populations have artificially sharpened peaks in this diagnostic due to their being better described by single exponentials. Notably also, in Figure 6, at low $\mathrm{[Fe/H]}$ the inner profiles flatten faster than the outer profile, whereas the higher $\mathrm{[Fe/H]}$ populations show the opposite behaviour. For example, in the $-0.3<\mathrm{[Fe/H]}<-0.2\ \mathrm{dex}$ bin, the outer profile appears to remain roughly constant in slope between 1 and 6 Gyr, while the inner profile flattens significantly. The opposite is seen in the $0.1<\mathrm{[Fe/H]}<0.2\ \mathrm{dex}$ bin, where the outer profile flattens considerably with age.

4.2 The vertical profile of the disk

We now examine the variation of $h_{Z}$ as a function of radius in mono-age, mono-$\mathrm{[Fe/H]}$ populations. Because mono-age, mono-$\mathrm{[Fe/H]}$ populations are well described by a single scale height, which is modified by a flaring term $R_{\text{flare}}$, this means that $h_{Z}$ is weakly dependent on $R$ for profiles which flare. We show vertical profiles for age-$\mathrm{[Fe/H]}$ bins with $>30$ stars in Figure 7, adopting the same shading as Figure 6 to draw the eye to the profiles with greater mass contribution, and adding a dashed line representing $0.3$ kpc for reference.

Figure 7 suggests that the disk is thicker as traced by older populations. All $\mathrm{[Fe/H]}$ bins show a thickening as age increases. This is clear in the left panel of Figure 9, which shows the surface-mass density weighted mean variation of $h_{Z}$ with age. The mean $h_{Z}$ spans the range between 0.8 and 0.2 kpc. The high $\mathrm{[}\alpha\mathrm{/Fe]}$ populations show a bump in the mean $h_{Z}$ at 8 Gyr, but $h_{Z}$ generally increases with age, similarly to the low $\mathrm{[}\alpha\mathrm{/Fe]}$ populations. The shapes of the profiles of the youngest populations in Figure 7 in the low $\mathrm{[}\alpha\mathrm{/Fe]}$ subsample show little variation with $\mathrm{[Fe/H]}$, and this trend generally continues to older ages. This is also reflected in the low uncertainties associated with the blue points in the left panel of Figure 9.

The high $\mathrm{[}\alpha\mathrm{/Fe]}$ profiles are generally flat, indicating that these populations show little flaring. By multiplying together the PDFs of the posterior distribution of fits to $R_{\mathrm{flare}}$, we determine that the high $\mathrm{[}\alpha\mathrm{/Fe]}$ populations have an average $R_{\mathrm{flare}}^{-1}=-0.06\pm 0.02$. The low alpha populations flare more strongly, with an average $R_{\mathrm{flare}}^{-1}=-0.12\pm 0.01$. There is, however, some variation in the flaring as a function of age, so it may not be sensible to ascribe a single $R_{\mathrm{flare}}^{-1}$ to all the populations. We show the variation of $R_{\mathrm{flare}}^{-1}$ with age in the right panel of Figure 9, showing $R_{\mathrm{flare}}^{-1}$ rather than $R_{\mathrm{flare}}$ so that values close to 0 are represented properly. The surface-mass density weighted mean $R_{\mathrm{flare}}^{-1}$ of the low $\mathrm{[}\alpha\mathrm{/Fe]}$ populations increases as a function of age, meaning that the most flared populations are the youngest. The behaviour appears opposite for the high $\mathrm{[}\alpha\mathrm{/Fe]}$ populations, whose mean $R_{\mathrm{flare}}^{-1}$ seems to decrease with age, but this is determined with low significance as $R_{\mathrm{flare}}^{-1}$ measurements are noisier for these populations.

4.3 The mass contribution of mono-age, mono-$\mathrm{[Fe/H]}$ populations

We now present the results from the calculation of the surface-mass density at the solar radius using the method described in Section 3.3. We compute surface-mass density $\Sigma_{R_{0}}$ estimates for each age-$\mathrm{[Fe/H]}$ bin, for the high and low $\mathrm{[}\alpha\mathrm{/Fe]}$ samples. When quoting the surface-mass densities, we also quote estimates of the systematic uncertainties. We evaluate the sources of these uncertainties in Section 5.1.

We combine the mass contributions of the high and low $\mathrm{[}\alpha\mathrm{/Fe]}$ mono-age, mono-$\mathrm{[Fe/H]}$ populations, and plot the estimates as a function of age and $\mathrm{[Fe/H]}$ in Figure 10. This figure essentially represents the mass-weighted age-$\mathrm{[Fe/H]}$ distribution at the solar radius, that is, the probability distribution for age and $\mathrm{[Fe/H]}$ for a randomly selected mass element. The distribution varies smoothly with no sharp peaks, and the surface-mass density increases linearly with both age and $\mathrm{[Fe/H]}$, peaking at $1<\mathrm{age}<3$ Gyr, $0.0<\mathrm{[Fe/H]}<0.1$ dex. The mass increases more smoothly with $\mathrm{[Fe/H]}$ than with age, but there is little mass in the highest $\mathrm{[Fe/H]}$ bin, creating a ridge in the marginalised distribution. The marginalised distributions as a function of age and $\mathrm{[Fe/H]}$ show no sign of bimodality, and there is little sign of a bimodality in age at fixed $\mathrm{[Fe/H]}$. It should be mentioned again here that the age uncertainties may be larger than the bin width, particularly in older bins, which would cause an artificial blurring of a density edge in teh distribution along the age axis. Therefore, we cannot presently determine to high significance that there are no discontinuities in this distribution.

An alternative way to look at the surface-mass density distributions is to retain the division in $\mathrm{[}\alpha\mathrm{/Fe]}$. We find that the high $\mathrm{[}\alpha\mathrm{/Fe]}$ populations contribute $\Sigma_{R_{0},\ \mathrm{tot}}=3.0_{-0.5}^{+0.4}\mathrm{(stat.)}_{-0.6}^{+0.6}\mathrm{(syst.)}\ \mathrm{M_{\odot}\ pc^{-2}}$ to the total surface-mass density at the solar radius, whereas the low $\mathrm{[}\alpha\mathrm{/Fe]}$ populations contribute $\Sigma_{R_{0},\ \mathrm{tot}}=17.1_{-2.4}^{+2.0}\mathrm{(stat.)}_{-1.9}^{+4.4}\mathrm{(syst.)}\ \mathrm{M_{\odot}\ pc^{-2}}$, giving a total surface-mass density in stars at $R_{0}$ of $\Sigma_{R_{0},\ \mathrm{tot}}=20.0_{-2.9}^{+2.4}\mathrm{(stat.)}_{-2.4}^{+5.0}\mathrm{(syst.)}\ \mathrm{M_{\odot}\ pc^{-2}}$. We plot the individual surface-mass contributions for the separated low and high $\mathrm{[}\alpha\mathrm{/Fe]}$ populations in Figure 11, adopting different color scales in each panel, to highlight the behaviour of the high $\mathrm{[}\alpha\mathrm{/Fe]}$ populations, which contribute little mass in comparison to the low $\mathrm{[}\alpha\mathrm{/Fe]}$. The low $\mathrm{[}\alpha\mathrm{/Fe]}$ mass is mostly concentrated at young age and towards higher $\mathrm{[Fe/H]}$, although there is mass even at the oldest ages. The high $\mathrm{[}\alpha\mathrm{/Fe]}$ mass is concentrated towards older ages, but interestingly the distribution extends to high $\mathrm{[Fe/H]}$, and we detect mass at some $\mathrm{[Fe/H]}$ in every age bin, but at much lower levels. The tails of the distributions of the low and high $\mathrm{[}\alpha\mathrm{/Fe]}$ populations overlap somewhat in age-$\mathrm{[Fe/H]}$ space, around 6 Gyr ago, and there is a hint of a sequence extending from old, low $\mathrm{[Fe/H]}$ and high $\mathrm{[}\alpha\mathrm{/Fe]}$ populations, to young, high $\mathrm{[Fe/H]}$ and low $\mathrm{[}\alpha\mathrm{/Fe]}$ populations, which is somewhat visible in the combined histogram. There is no clear bimodality in age at fixed $\mathrm{[Fe/H]}$ in the combined histogram, owing to the very low mass contribution of the old, high $\mathrm{[}\alpha\mathrm{/Fe]}$ populations.

We have established that the vertical spatial distributions of mono-age, mono-$\mathrm{[Fe/H]}$ populations are well described by single exponentials with characteristic $h_{Z}$. Next, we use this information to generate the mass-weighted distribution of $h_{Z}$, which is representative of the probability distribution function for $h_{Z}$, $p(h_{Z})$. For a random stellar mass element, this function gives the probability density for the $h_{Z}$ of the component to which it belongs. We show this relation in Figure 12, where coloured points represent the individual density contributions of mono-age, mono-$\mathrm{[Fe/H]}$ populations, and the coloured histograms their co-addition within $\sim 0.1$ kpc wide bins in $h_{Z}$ for the low and high $\mathrm{[}\alpha\mathrm{/Fe]}$ populations (purple and green, respectively). The dashed histogram represents the resulting total $p(h_{Z})$. Scatter points are coloured by the age of the population they represent. The total distribution is smooth, resulting from the superposition of the low and high $\mathrm{[}\alpha\mathrm{/Fe]}$ distributions, which overlap significantly. The total $\Sigma_{R_{0}}$ (dashed histogram) declines exponentially with $h_{Z}$, and is unimodal with no gaps. The trends of both $h_{Z}$ and $\Sigma_{R_{0}}$ with age seen in Figures 8 and 10 are recovered here, although it is surprising that the trend of $h_{Z}$ with age at the high $h_{Z}$ end does not appear as obvious here.

We can also now mass-weight and combine the fitted density profiles to attain the surface-mass density profile of the Milky Way as a function of age, $\mathrm{[Fe/H]}$ and $\mathrm{[}\alpha\mathrm{/Fe]}$. The resulting profiles are displayed in Figure 13. The different nature of the low and high $\mathrm{[}\alpha\mathrm{/Fe]}$ populations in terms of spatial structure is clear here, with the low $\mathrm{[}\alpha\mathrm{/Fe]}$ profile having a clear break between 8 and 10 kpc, and the high $\mathrm{[}\alpha\mathrm{/Fe]}$ declining exponentially with $R$. It is interesting to note that extrapolation by eye of the high and low $\mathrm{[}\alpha\mathrm{/Fe]}$ profiles to low $R$ would result in the high $\mathrm{[}\alpha\mathrm{/Fe]}$ population becoming dominant over the low. The total profile appears roughly flat out to $\sim 10$ kpc. However, we strongly emphasize that this is not determined to high significance, as even when only the uncertainties from the fitting procedure are included, one could describe the profile as exponentially declining with $R$ within $R<R_{0}$. The inclusion of the other sources of uncertainty on the surface-mass density estimates would further decrease the significance of the apparent flattening. For example, the systematic uncertainties (discussed in Section 5.1) act to increase the fraction of surface-mass density contributed by the high $\mathrm{[}\alpha\mathrm{/Fe]}$ populations, which would only increase the slope of the inner exponential. Using dynamcial tracers, Bovy & Rix (2013) find that the surface density should decline exponentially with R, so it seems logical to assume that the inner profile should not be increasing with R.

As a function of age, the peak in the surface-mass density visible in the youngest population becomes less prominent, and the profile becomes a roughly single exponential at the oldest ages (i.e., it monotonically decreases with $R$). The behaviour with $\mathrm{[Fe/H]}$ is more complex, but the variation of the peak radius with $\mathrm{[Fe/H]}$ is obvious, and the turnover in the total profile at $\sim 10$ kpc appears to be a result of the outermost breaks.

5 Discussion

In the above analysis, we have, for the first time, determined the detailed structure of the Milky Way’s disk as a function of stellar age and $\mathrm{[Fe/H]}$. In our method, we have drawn heavily from previous dissections of the disk into its mono-abundance constituents (MAPs; Bovy et al., 2012b; Bovy et al., 2016b), and so use these previous findings as a benchmark with which to compare these results. We also show that our results are also broadly consistent with other measurements, whilst shedding new light onto the problem of the formation of the Milky Way disk.

5.1 Surface-mass density systematics

We first adress the sources of systematic uncertainty in our surface-mass density estimates, which are pertinent to the following discussions. Our total local surface-mass density, including the correction for stars missing from the age catalogue, but before accounting for any other systematic uncertainties, is $\Sigma_{R_{0},\text{tot}}=20.0_{-2.9}^{+2.4}\ \mathrm{M_{\odot}\ pc^{-2}}$. This result is roughly two thirds as large as previous estimates, which are of order $\Sigma_{R_{0},\text{tot}}\sim 30\ \mathrm{M_{\odot}\ pc^{-2}}$ (e.g. Flynn et al., 2006; Bovy et al., 2012a; McKee et al., 2015). From canonical stellar evolution it is known that giants contribute very little to the total stellar mass in any population. For instance, McKee et al. (2015) find that giants make up $\sim 2\%$ of the local stellar mass. By virtue of this fact, conversions of the stellar mass inferred from giant-star counts to that of the total underlying stellar population require a multiplication of the observed counts by a factor of $\sim 50$, meaning that any uncertainty in the star counts is amplified in the final surface-mass density estimate. Our quoted statistical error estimates, however, which account for Poisson fluctuations in the stellar counts, cannot fully account for the discrepancy.

We first evaluate whether such a discrepancy may be due to the assumed IMF or stellar evolution model. Tests adopting exponential Chabrier (2003) and Kroupa (2001) IMFs for the mass calculation resulted in variations of the final $\Sigma_{R_{0},\text{tot}}$ estimate of the order $\sim 1\ \mathrm{M_{\odot}\ pc^{-2}}$, which we incorporate into the systematic error budget. We also re-ran our analysis on the basis of the BaSTI stellar evolution models (Pietrinferni et al., 2004), for which there also exists calculations for $\alpha$-enhanced stars (Pietrinferni et al., 2006), which produced comparable estimates to the PARSEC models (after correcting for the fact that the lowest mass in the BaSTI isochrones is $0.5M_{\odot}$ as opposed to $0.1M_{\odot}$ in the PARSEC models). We also compute the mass using only APOGEE fields away from the plane (with $|b|\geq 6^{\circ}$), to test for the effects of extinction on the star counts, but attain results within the Poisson uncertainties of the original estimate.

We also apply our analysis procedure to a basic Monte Carlo mock sample, to check the method for converting observed counts to the real number density $N(R_{0})$. We sample stars on a broken exponential density distribution with exponential flare then select points within APOGEE fields out to an imposed distance cut (which allows a simple reconstruction of the selection, and calculation of the effective volume). We calculate $N(R_{0})$ analytically, and via our method, and find results which are consistent with the input parameters of the model broken exponential profiles, within the Poisson errors, for a wide variety of input parameters.

As mentioned in Section 2.3, we find that, after making corrections to the ages, that the model returns ages greater than 13 Gyr for a sizeable number of stars (Martig et al., 2016b, limit ages to 13 Gyr in their table). While these stars make up approximately $10\%$ of the final sample (3020 stars), they are not included in the number counts in each mono-age mono-$\mathrm{[Fe/H]}$ bin for calculation of the surface-mass density. Adding an extra $10\%$ of counts to each bin (in the same way as the extra $25\%$ is added in Section 2.3) introduces an extra systematic uncertainty of roughly $1\mathrm{M_{\odot}\ pc^{-2}}$ in each $\mathrm{[}\alpha\mathrm{/Fe]}$ sub-sample. However, readers should take into account that this simple correction does not account for a scenario where the stars with ages fitted $>13$ Gyr might have a specific distribution in age, casting more counts in some bins (which might have more mass contribution per star) than others. For example, if these stars were all old, then the actual surface mass density in older bins would be higher than that found here, which would increase the total surface-mass density estimate.

From the above, we conclude that the majority of the systematic discrepancy is likely not due to the assumed IMF, stellar evolution model, dust extinction, or some peculiarity in the age measurements which affects star counts in the bins used. At this stage, it is difficult to understand what is the possible origin of this discrepancy with other works in the literature. Interestingly, our study is the only one employing giant stars as the stellar population tracer, which may point to possible systematics in the theoretical isochrones, or the APOGEE stellar parameters, or a combination thereof. It has recently been demonstrated by Masseron & Hawkins (2017) that there may be significant issues with the spectroscopic determination of stellar surface gravity, which is dependent on the star’s evolutionary state. We find some discrepancy between the $\log{g}$ of the red clump between the PARSEC isochrones and the data, of the order $\sim 0.2$ to $0.3$ dex (similar to that found by Masseron & Hawkins, 2017, albeit based on APOGEE-DR13 data), which could conceivably lead to problems in our conversion. We test for the effect of systematics in the $\log{g}$ scale, shifting the $\log{g}$ cut of the isochrones to lower and higher $\log{g}$ by $0.3$ dex. We find that shifting the $\log{g}$ cut by $-0.3$ dex increases the surface-mass density estimate by $5.0\ \mathrm{M_{\odot}\ pc^{-2}}$. Increasing the $\log{g}$ cut by $0.3$dex results in a decrease of $2.4\ \mathrm{M_{\odot}\ pc^{-2}}$. It therefore seems plausible that the discrepancy results from a systematic difference between the $\log{g}$ scales of the theoretical isochrones and APOGEE, and so we incorporate these shifts into the systematic error estimate.

Upon inclusion of the systematic uncertainties from IMF variations and differences in the surface gravity scales, we attain a final estimate of the total local surface-mass density in visible stars of $\Sigma_{R_{0},\text{tot}}=20.0_{-2.9}^{+2.4}\mathrm{(stat.)}_{-2.4}^{+5.0}\mathrm{(syst.)}\ \mathrm{M_{\odot}\ pc^{-2}}$, from the addition of the low $\mathrm{[}\alpha\mathrm{/Fe]}$ surface-mass density of $\Sigma_{R_{0},\text{tot}}=17.1_{-2.4}^{+2.0}\mathrm{(stat.)}_{-1.9}^{+4.4}\mathrm{(syst.)}\ \mathrm{M_{\odot}\ pc^{-2}}$, and the high $\mathrm{[}\alpha\mathrm{/Fe]}$ value of $\Sigma_{R_{0},\text{tot}}=3.0_{-0.5}^{+0.4}\mathrm{(stat.)}_{-0.6}^{+0.6}\mathrm{(syst.)}\ \mathrm{M_{\odot}\ pc^{-2}}$. If the $\log{g}$ systematics are as large as $-0.3$ dex, then our result is in agreement with the recent estimate from McKee et al. (2015), of $27\pm 2.7\ \mathrm{M_{\odot}\ pc^{-2}}$.

A recent compilation of measurements of the thick and thin disks found that the thick-thin disk surface density ratio at $R_{0}$ is $f_{\Sigma}=15\%\pm 6\%$ (Bland-Hawthorn & Gerhard, 2016). Our results find $f_{\Sigma}=18\%\pm 5\%$ for high-low $\mathrm{[}\alpha\mathrm{/Fe]}$ disk surface-mass density ratio, consistent with that estimate. While a better understanding of the possible systematics between the theoretical isochrones and APOGEE is beyond the scope of this paper, we have shown that even a slight difference in the $\log{g}$ scale can bring our results in line with existing estimates. This suggests that our surface-mass density measurement discrepancy is indeed systematic, and that the high and low $\mathrm{[}\alpha\mathrm{/Fe]}$ disks may still have some relation to the thick and thin components measured by these studies, which are mainly based on geometric decompositions of the disk.

5.2 Comparison with MAPs results

We first discuss our density fits in comparison to the MAP measurements of Bovy et al. (2012b); Bovy et al. (2016b). Such a comparison is important because our method is based on an extension of that developed by Bovy et al. (2016b) to the case of RGB stars, whose distances are far more uncertain than those of RC stars. Bovy et al. (2016b) used the APOGEE RC sample, to find the structure of populations in narrow bins of $\mathrm{[}\alpha\mathrm{/Fe]}$ and $\mathrm{[Fe/H]}$, or MAPs. These MAPs represent stellar populations with a distribution of ages, but their interpretation assumes a significant relationship between age, $\mathrm{[}\alpha\mathrm{/Fe]}$ and $\mathrm{[Fe/H]}$. We can now compare results when the third parameter, age, is known.

Bovy et al. (2016b) showed that the radial distribution of low $\mathrm{[}\alpha\mathrm{/Fe]}$ MAPs is well described by a broken exponential, and we confirm this result, showing that each of the low $\mathrm{[}\alpha\mathrm{/Fe]}$, mono-age, mono-metallicity populations is also described by a radially broken exponential. We also show that the older, high $\mathrm{[}\alpha\mathrm{/Fe]}$ populations are instead described by a single exponential, which is in good agreement with the findings of Bovy et al. (2016b).

The dependence of the radial distribution of mono-$\mathrm{[Fe/H]}$ populations on age is interesting in this regard. The low $\mathrm{[}\alpha\mathrm{/Fe]}$ population, for which our sample covers a wide range of ages with high signal-to-noise, shows a broadening of the profile around a density peak towards older populations, at all $\mathrm{[Fe/H]}$. This effect does not appear to be present in the high $\mathrm{[}\alpha\mathrm{/Fe]}$ population (although some populations have slight evidence of a break at low significance), which suggests that it was formed and evolved differently. We discuss the implications of these findings in Section 5.4.

We also confirm the results of Bovy et al. (2016b) which showed that the break radius, $R_{\mathrm{peak}}$ is a declining function of $\mathrm{[Fe/H]}$. We show that, in the low $\mathrm{[}\alpha\mathrm{/Fe]}$ population, at fixed age, $R_{\mathrm{peak}}$ moves to smaller radii as $\mathrm{[Fe/H]}$ increases. As a function of age, the amplitude of this variation increases. The difference in $R_{\mathrm{peak}}$ at the highest and lowest $\mathrm{[Fe/H]}$ in the 6 Gyr bin is $\sim 6$ kpc, which is identical with that of the profiles shown in Figure 11 of Bovy et al. (2016b).

Bovy et al. (2016b) found that low $\mathrm{[}\alpha\mathrm{/Fe]}$ MAPs were fit well with a $h_{Z}(R)$ which was slowly exponentially flaring with $R_{\mathrm{flare}}^{-1}=-0.12\pm 0.01\ \mathrm{kpc^{-1}}$. We confirm this result, finding also that low $\mathrm{[}\alpha\mathrm{/Fe]}$ populations have, on average, $R_{\mathrm{flare}}^{-1}=-0.12\pm 0.01\ \mathrm{kpc^{-1}}$, but we also find that $R_{\mathrm{flare}}^{-1}$ shows considerable variation with age (around this mean) in the low $\mathrm{[}\alpha\mathrm{/Fe]}$ populations. While Bovy et al. (2016b) found that high $\mathrm{[}\alpha\mathrm{/Fe]}$ populations were consistent with a having $R_{\mathrm{flare}}^{-1}=0.0\pm 0.02\ \mathrm{kpc^{-1}}$, we find that these populations have $R_{\mathrm{flare}}^{-1}=-0.06\pm 0.02\ \mathrm{kpc^{-1}}$, showing some evidence of flaring, albeit at a lower level and lower significance than the low $\mathrm{[}\alpha\mathrm{/Fe]}$ populations.

It was also shown by Bovy et al. (2012b); Bovy et al. (2016b) that the $h_{Z}$($\mathrm{[}\alpha\mathrm{/Fe]}$,$\mathrm{[Fe/H]}$) of MAPs smoothly spans the range between 0.2 and 1 kpc. We also confirm and extend this result, showing that $h_{Z}$(age,$\mathrm{[}\alpha\mathrm{/Fe]}$,$\mathrm{[Fe/H]}$) varies smoothly between a maximum $h_{Z}$ of $\sim 1.2$ kpc in the high $\mathrm{[}\alpha\mathrm{/Fe]}$, low $\mathrm{[Fe/H]}$, older populations, down to a minimum of $\sim 0.2$ kpc in the youngest, low $\mathrm{[}\alpha\mathrm{/Fe]}$, $\mathrm{[Fe/H]}$ rich populations. We can also directly compare our Figure 12 with Figure 2 of Bovy et al. (2012a), which showed that the mass-weighted $h_{Z}$ distribution is not bimodal but smoothly declines with $h_{Z}$. We confirm that result, showing that for mono-age, mono-metallicity populations, the mass weighted $h_{Z}$ distribution shows no sign of bimodality; The low and high $\mathrm{[}\alpha\mathrm{/Fe]}$ populations’ $h_{Z}$ distributions are distinct but overlap significantly, generating a smooth distribution. This presents an interesting new look at the interplay between spatial and chemical structure in the disk, as there is a clear mixing spatially of the two chemically separated populations. The implications of this finding are intriguing, and we discuss them further in Section 5.4.

Bovy et al. (2012a) also made a measurement of the local surface-mass density, finding (in the original application of the method used here) that SEGUE G-type dwarfs yield an estimate of $\Sigma_{R_{0},\text{tot}}=30\pm 1\ \mathrm{M_{\odot}\ pc^{-2}}$, which is in good agreement with other studies based on different samples (e.g. Flynn et al., 2006; McKee et al., 2015). Comparatively, our result is somewhat smaller, even when systematics are taken into account. There are a number of differences between this study and that of Bovy et al. (2012a). For example, the increased radial coverage, adoption of RGB stars as a tracer, and the fits based on mono-age, mono-metallicity (rather than mono-abundance) populations. However, as mentioned above in Section 5.1, our results, after accounting for systematic uncertainties, appear in good agreement with other, more recent estimates (McKee et al., 2015).

5.3 Comparison with other Milky Way disk studies

We now compare qualitatively the findings of our analysis with the broader body of knowledge regarding the Milky Way disk structure (see, e.g. Rix & Bovy, 2013; Bland-Hawthorn & Gerhard, 2016, for recent reviews) . In comparing our results with those from previous studies, we are constrained to making mostly qualitative considerations, as previous work is based on fits of single exponentials to the radial component of the stellar density distribution.

This work strongly constrains the structure of both the low and high $\mathrm{[}\alpha\mathrm{/Fe]}$ components in the Galactic disk, which are commonly considered to be interchangeable with the thin and thick components (as asserted by, e.g. Bensby et al., 2004; Adibekyan et al., 2012; Fuhrmann, 1998). We have shown that the $\mathrm{[}\alpha\mathrm{/Fe]}$ rich component, while corresponding to a thicker configuration in general, is the product of individual mono-age populations of varying thickness. We find that the $\mathrm{[}\alpha\mathrm{/Fe]}$ rich populations span the range $0.4<h_{Z}<1$ kpc, with $h_{Z}$ increasing with age. Studies of the vertical disk structure which fit a double exponential find a thick disk scale height of $\sim 1$ kpc (e.g. Gilmore & Reid, 1983; Jurić et al., 2008), which is fully consistent with measurements of the thickest high $\mathrm{[}\alpha\mathrm{/Fe]}$, old, mono-age populations in our analysis. However, we again stress here that the age uncertainties at old ages may be significantly larger than the bin size, which may cause a blurring of these trends, and should be accounted for when comparing these results to models. As an example, in the worst case scenario, assuming gaussian errors, the oldest bins (between 7 and 13 Gyr) may be contaminated by up to $50\%$ of the stars which should be assigned to neighbouring bins, at the oldest end, with the fraction dropping off quickly at younger ages. We briefly discuss the implications of the worst case blurring on our interpretation of these trends in Appendix B

Regarding the radial scale length of the thick component, we find an obvious discrepancy with literature values, whereby our thick, high $\mathrm{[}\alpha\mathrm{/Fe]}$, populations have an average $h_{R,\text{out}}=1.9\pm 0.1$ kpc, while the aforementioned studies, who define the thick disk geometrically, find values of the order $\sim 4$ kpc (Ojha, 2001; Jurić et al., 2008). This discrepancy appears to arise in the choice of definition of the measured population between a geometric or chemical abundance selection, with many studies finding a scale length for the abundance-selected $\alpha$-rich disk in the range $h_{R}=2.0\pm 0.2$ kpc (Bland-Hawthorn & Gerhard, 2016; Bovy et al., 2016b; Bovy et al., 2012b; Cheng et al., 2012b). It should be noted here that our $h_{R,\text{out}}$ would likely be in even better agreement with this value, had we accounted for the 5% systematic discrepancy in the distances (shown in Figure 3). Martig et al. (2016a) recently showed evidence for a radial age gradient in the Milky Way, suggesting this as a source of disagreement between abundance-selected and geometric studies of the thick disk components, where the geometric selected studies see an extended thick disk which is made up of flared low $\mathrm{[}\alpha\mathrm{/Fe]}$ populations.

We also examine claims of a sharp decline in the stellar density at $R\sim 13.5$ kpc (e.g. Reylé et al., 2009; Sale et al., 2010) in light of our results. Sale et al. (2010) fit a single exponential density profile with scale-length $\sim 3$ kpc to A stars (which preferentially selects stars younger than $\sim 100$ Myr old) and found that after $R\sim 13$ kpc, a model with shorter scale-length was necessary to explain the increased rate of decline in stellar density. Our total profile in Figure 13 begins to decline after $R\sim 10$ kpc. The uncertainties in the measurement of the mass contribution of each profile may cause some discrepancy here, as implying a higher mass on older or more metal poor populations would shift this turn-over to higher radii. We also fit older populations than Sale et al. (2010), which, under the inside-out formation paradigm, might suggest another reason for such a discrepancy, as older populations would be more centrally concentrated. We confirm the assertion of Bovy et al. (2016b) that this break, clearly visible in the total stellar distribution, is attributable to the outermost break of the mono-$\mathrm{[Fe/H]}$ profiles - which are shown in the bottom panel of Figure 13. External disk galaxies are also observed to have such a truncation in their stellar density profiles (e.g. Pohlen & Trujillo, 2006).

5.4 Implications for the formation of the Galactic disk

In light of the above discussion, we now present the implications of our results for the formation of the Galactic disk. In this paper, we present a detailed dissection of the disk by age, $\mathrm{[Fe/H]}$,and $\mathrm{[}\alpha\mathrm{/Fe]}$, and as such, present a previously unseen picture of the dominant structure of the Milky Way. In studying mono-age populations, we can perform a more direct comparison than previously possible with numerical simulations of Milky Way type galaxies, which tend to use age information in the absence of detailed chemical modelling.

5.4.1 Disk flaring, profile broadening and radial migration

By estimating the density profiles of mono-age populations, we place novel constraints on radial migration and its effects on the structure and evolution of the disk. We have two key observables which provide this insight: the flaring of the disk, which has been considered as an effect of vertical action conserving radial migration (where stars have greater vertical excursions as they migrate outward, e.g. Minchev et al., 2012), and the broadening of the density profiles around the peak with time, which we discuss as a potential new indicator of radial migration. The right panel of Figure 9 shows a clear trend of increasing $R_{\mathrm{flare}}^{-1}$ with age such that the youngest populations flare most. This behaviour is distinct from the results of Bovy et al. (2016b), which found that low $\mathrm{[}\alpha\mathrm{/Fe]}$ populations were described by a single $R_{\mathrm{flare}}^{-1}$. It is, however, conceivable that if low $\mathrm{[}\alpha\mathrm{/Fe]}$ populations of all ages are combined, the resulting population may have a similar behaviour as that in Bovy et al. (2016b). Indeed, we find an average $R_{\mathrm{flare}}^{-1}=-0.12\pm 0.01\ \mathrm{kpc^{-1}}$, which is in good agreement with the value from Bovy et al. (2016b). If populations become more flared as radial migration proceeds, then it becomes difficult to reconcile our result with that of a disk whose stars continually underwent radial migration, largely unperturbed by any mergers that might cause structural discontinuity (e.g. Martig et al., 2014a), especially under suggestions that mergers actually reduce flaring from radial migration (e.g. Minchev et al., 2014a). In this context, this result is indicative of an old population in the Milky Way which has undergone some mergers, reducing the flaring in the oldest populations. It should be noted here, however, that the age uncertainties (which can be as large as $40\%$) could effectively artificially increase the age bin size, super-imposing populations with different scale heights and flare, and reducing the overall flaring profile.

We showed in Figures 8 and 6 that the radial surface density profiles of low $\mathrm{[}\alpha\mathrm{/Fe]}$ populations become smoother with age. Interestingly, the position of the break radius does not vary monotonically with age. Assuming that the peak radius is at the equilibrium point of chemical evolution for a given population (where the consumption of gas and its dilution are balanced, as discussed in Bovy et al., 2016b), then one might consider that such a broadening would occur if stars that formed near the equilibrium point migrated inwards and outwards over time. If these assumptions are correct, then the specific surface density profile shapes might provide insights into how radial migration has proceeded in the disk. For example, if the slope of the inner or outer profile change slope differently, this might suggest that migration has been asymmetric (i.e. more mass has moved in than out or vice versa). Comparing the $-0.2<\mathrm{[Fe/H]}<-0.1$ dex and $0.0<\mathrm{[Fe/H]}<0.1$ dex bins in Figure 6, it seems that the inner slope of the former profile decreases more with age, whereas the outer slope of the latter decreases more strongly, suggesting that the former population has preferentially migrated in, whereas the latter migrated out. it is important to point out that under the interpretation that the increasing profile width is due to migration efficiency, we make an assumption that stars of a given $\mathrm{[Fe/H]}$ and $\mathrm{[}\alpha\mathrm{/Fe]}$ must have been born with the same width throughout cosmic time (as discussed by, e.g. Minchev et al., 2017).

This picture is also consistent with the suggestion by, e.g, Hayden et al. (2015); Loebman et al. (2011), that the changing skew in the MDF as a function of Galactocentric radius is caused by such a mechanism. We show the mass-weighted $\mathrm{[Fe/H]}$ distribution for the low $\mathrm{[}\alpha\mathrm{/Fe]}$ fits at different Galactocentric radii in Figure 14, showing that our results both find a radial metallicity gradient, and qualitatively reproduce the skew found by Hayden et al. (2015). Unlike Hayden et al. (2015), our analysis fully corrects for sample-selection and stellar-population biases in reconstructing the MDF. Grand et al. (2016) also found similar behaviour in a simulated galaxy, finding that spiral structure induces different migration patterns, dependent on birth radius.

Interestingly, we detect the possible flaring of old, high $\mathrm{[}\alpha\mathrm{/Fe]}$ populations, with average $R_{\mathrm{flare}}^{-1}=-0.06\pm 0.02\ \mathrm{kpc^{-1}}$. However, the right panel of Figure 9 shows that the flaring of the high $\mathrm{[}\alpha\mathrm{/Fe]}$ populations does not vary as strongly with age as the low $\mathrm{[}\alpha\mathrm{/Fe]}$ populations (although it should be noted that most of the mass in the high $\mathrm{[}\alpha\mathrm{/Fe]}$ populations is concentrated at older age anyway). The detection of even a slight flare in these populations is surprising, as Bovy et al. (2016b) found that the high $\mathrm{[}\alpha\mathrm{/Fe]}$ MAPs did not have flare. Again, this may be an effect of the superposition of multiple mono-age populations within the MAPs. Minchev et al. (2015), for example, found that co-eval populations in simulated galaxies always flare, and suggested that the superposition of such flares might be an explanation for thickened disk components. Minchev et al. (2017) showed that the superposition of mono-age populations within MAPs can introduce decreased flaring in high $\mathrm{[}\alpha\mathrm{/Fe]}$ populations, whilst the mono-age populations themselves still flare. Comparison of these results with Bovy et al. (2016b) seems to present a consistent scenario. It should be noted here however, that Stinson et al. (2013) found that MAPs in their simulation were coeval in general.

Our results show that the oldest populations are thicker, centrally concentrated, and display the least flaring, whilst the youngest populations, which show the most flaring, have the thinnest vertical distribution (smallest $h_{Z}$). Between these extremes, consecutive populations in age form a continuum, when the combined low and high $\mathrm{[}\alpha\mathrm{/Fe]}$ structure is considered (see Figure 9). It is clearly conceivable then, that the geometrically defined thick disk, found to have large scale-length (e.g., Jayaraman et al., 2013; Jurić et al., 2008), may be the superposition of these flared (young) and naturally thick (old) components. An obvious consequence of this scenario would be an age-gradient at high $Z$ above the disk plane, which has been recently shown to be present in the APOGEE data by Martig et al. (2016a). It was also recently shown that in the Gaia-ESO survey data, the mean structural characteristics of the abundance selected thick and thin disks appear to overlap at $\mathrm{[M/H]}\sim-0.25$ dex and $\mathrm{[}\alpha\mathrm{/Fe]}$ $\sim 0.1$ dex (Recio-Blanco et al., 2014), further presenting a scenario where the thick and thin disk components are not necessarily separable from one another, or at least not in abundance space.

5.4.2 Inside out formation, and the overall vertical disk structure

The formation of the Galactic disk is commonly framed in the paradigm of inside-out formation (e.g., Bird et al., 2013; Kobayashi & Nakasato, 2011; Larson, 1976; Matteucci & Francois, 1989). More recently, the effects of radial migration (e.g. Sellwood & Binney, 2002) were added, in order to produce models that agree better with the observations (e.g. Loebman et al., 2011; Schönrich & Binney, 2009a; Spitoni et al., 2015; Kubryk et al., 2015) . Our measurements of the peak radius of mono-age populations place strong empirical constraints on the evolution of the Milky Way disk over time. The behaviour of $R_{\mathrm{peak}}$ with age and $\mathrm{[Fe/H]}$ is shown in Figure 15. We find that the surface-mass density weighted mean $R_{\mathrm{peak}}$ of low $\mathrm{[}\alpha\mathrm{/Fe]}$ populations remains roughly constant with age, whilst the dispersion about the mean increases with age. This finding is qualitatively consistent with that of e.g. Anders et al. (2016b), who show that the radial metallicity gradient decreases with the age of the population considered. Our results show that the density peaks of mono-$\mathrm{[Fe/H]}$ populations become more separated with age. As the mean $\mathrm{[Fe/H]}$ at a given $R$ is dictated by the dominant population in stellar density at that radius, this indicates that we also see a shallowing gradient in $\mathrm{[Fe/H]}$ with age. This is reinforced by our finding (also shown in Figure 15) that the mean $R_{\mathrm{peak}}$ decreases with $\mathrm{[Fe/H]}$.

Our results also place strong constraints on models for the formation of the vertical disk structure. We have already discussed that the vertical structure of the disk is commonly framed as having two geometrically distinct (but overlapping) vertical components (e.g. Gilmore & Reid, 1983). Our results confirm previous work (e.g. Bovy et al., 2012a) that shows that this picture, while providing an acceptable description of the data when analysing the whole population, is not complete when individual populations, either abundance- or age-selected, are considered. We have shown (in Figure 12) that for a random mass element, the probability that it belongs to a population of a given $h_{Z}$ exponentially declines as $h_{Z}$ increases. There are no apparent breaks in this relation at the resolution which we measure it, strongly suggesting that the spatial vertical disk structure is continuous. This finding, in stark contrast to the distinct discontinuities seen in the chemical structure of the entire disk (e.g Nidever et al., 2014; Hayden et al., 2015), presents an interesting conundrum for galaxy formation theory. How is it possible that the spatial structure of the disk be smooth and continuous, whilst the chemical structure portrays a clear discontinuity?

Theoretical studies have, thus far, presented some clues as to how galaxy disks such as that of the Milky Way might form. As most studies fit single exponential radial profiles to simulated disks, quantitative comparisons are difficult, yet qualitative considerations can be made. Stinson et al. (2013) used a maximum likelihood method similar to Bovy et al. (2012b) to fit density profiles to MAPs in a simulated galaxy and found a continuous distribution of scale heights, but also found that their simulation showed a strongly geometrically distinct thick disk component. Bird et al. (2013) made detailed measurements of the mono-age populations in a high-resolution hydrodynamic Milky Way-like galaxy simulation, and found, similarly to our results, that their scale heights gradually decreased with time, while the scale lengths increased, with populations forming thick and retaining that thickness in an ’upside-down and inside-out’ disk formation. Our results also find evidence of flaring in the thick components (as discussed in Section 5.2), which may point towards some structural evolution (via a process such as radial migration) after their formation. However, work on simulations by Bournaud et al. (2009) suggests early, turbulent gas as the origin for thicker disk components. A flare in the gas disk of the Milky Way, associated with the stellar component, is also observed in numerous studies (e.g. Feast et al., 2014; Kalberla et al., 2014; Lozinskaya & Karadashev, 1963), suggesting a formation of the disk with structural parameters similar to its progenitor gas disk.

5.4.3 The age-$\mathrm{[Fe/H]}$ distribution and the evolution of the disk

We now discuss how the present day structural parameters, in combination with the emergent picture of the mass distribution in age-$\mathrm{[Fe/H]}$ space at the solar radius, might offer deeper insights into the formation and evolution of the disk.

We find, as may be expected, that the high $\mathrm{[}\alpha\mathrm{/Fe]}$ populations contribute the majority of their mass at the solar radius at ages older than $\sim 6$ Gyr, although the mass contribution by old stars is extremely low compared to the younger populations. The middle panel of Figure 13 shows that, if the populations follow the density models that we fit, then the older stars become more dominant closer to the Galactic centre, which is suggestive of a weak mean radial age gradient, in qualitative agreement with theoretical predictions (e.g. Minchev et al., 2015) and observations of the thick disk (e.g. Martig et al., 2016a). It is therefore not surprising that the bottom panel of Figure 13 shows a clear variation in mean metallicity, in agreement with findings in other works (e.g. Cheng et al., 2012a; Anders et al., 2016b; Hayden et al., 2015). Only with high resolution hydrodynamical simulations, which accurately reproduce the stellar populations in galaxies, will it be possible to reconstruct the right combination of star formation history and radial mixing that led to these age and metallicity gradients, to gain a better understanding of the details of their formation.

In Figure 11, an overlap in age-$\mathrm{[Fe/H]}$ space is visible between the high and low $\mathrm{[}\alpha\mathrm{/Fe]}$ populations. While there appears to be mass at many $\mathrm{[Fe/H]}$ bins in the old populations, the overlap in age occurs at intermediate $\mathrm{[Fe/H]}$ at the solar radius. Previous studies have found that the youngest stars in the high $\mathrm{[}\alpha\mathrm{/Fe]}$ sequence overlap in age with the oldest and most $\mathrm{[Fe/H]}$ poor stars in the low $\mathrm{[}\alpha\mathrm{/Fe]}$ population (e.g. Haywood et al., 2013), and our findings appear to be consistent with that result. It should however, be noted, that at least some of this overlap is likely caused by the age uncertainties, which can be as high as 40%. If the low $\mathrm{[}\alpha\mathrm{/Fe]}$ population emerged from the remnants of the high $\mathrm{[}\alpha\mathrm{/Fe]}$, then it is likely that some sort of infall event must have occurred to return the ISM to low $\mathrm{[Fe/H]}$ and low $\mathrm{[}\alpha\mathrm{/Fe]}$ before forming those stars (as expressed by, e.g., Chiappini et al., 1997). These scenarios are also discussed in the context of the APOGEE results by Nidever et al. (2014). To fully understand this, however, we will likely require a chemodynamical model which reproduces the bimodality in $\mathrm{[}\alpha\mathrm{/Fe]}$ at fixed $\mathrm{[Fe/H]}$.

6 Conclusions

We have performed the first detailed dissection of the stellar populations of the Milky Way disk in age, $\mathrm{[Fe/H]}$ and $\mathrm{[}\alpha\mathrm{/Fe]}$ space, bridging the gap between the detailed observational understanding of mono-abundance populations (e.g. Bovy et al., 2012b; Bovy et al., 2016b) and the plethora of studies of co-eval stellar populations in simulated galaxies (e.g Bird et al., 2013; Martig et al., 2014a; Stinson et al., 2013). We have placed novel constraints on models for the formation of the Milky Way disk by combining detailed density models fit to the mono-age, mono-$\mathrm{[Fe/H]}$ populations of the low and high $\mathrm{[}\alpha\mathrm{/Fe]}$ disk, with surface mass density contributions calculated on the basis of these density fits and stellar evolution models. We summarise our key results as follows:

• •

Radial and vertical profiles The mono-age, mono-$\mathrm{[Fe/H]}$ populations of the $\mathrm{[}\alpha\mathrm{/Fe]}$ poor disk are well fit by a radially broken exponential, with a peak radius, $R_{\mathrm{peak}}$, that varies as a function of age and $\mathrm{[Fe/H]}$. We find that the distance between $R_{\mathrm{peak}}$’s of the low and high $\mathrm{[Fe/H]}$ populations increases with age, which we interpret as evidence for a decreasing $\mathrm{[Fe/H]}$ gradient with time (e.g Anders et al., 2016b). The radial variation of the stellar surface density of the high $\mathrm{[}\alpha\mathrm{/Fe]}$ mono-age populations is found to have insignificant breaks, and they are better fit by a single exponential in this disk region. As these populations are the oldest, this may be a sign of the disk evolution washing out the density peak over time, or may point to a different formation scenario for high $\mathrm{[}\alpha\mathrm{/Fe]}$ stars, where no density peak ever existed. These findings are in good agreement with earlier studies of MAPs (Bovy et al., 2016b). We measure an average high $\mathrm{[}\alpha\mathrm{/Fe]}$ population scale length of $h_{R,\text{in}}=1.9\pm 0.1$ kpc, and find scale heights between 600 and 1000 pc, in good agreement with current measures of the $\mathrm{[}\alpha\mathrm{/Fe]}$ rich disk scale length and height (e.g. those outlined in Bland-Hawthorn & Gerhard, 2016).

• •

Profile Broadening We show that the radial surface density profile of the low $\mathrm{[}\alpha\mathrm{/Fe]}$ populations broadens with age in a given $\mathrm{[Fe/H]}$ bin, which we interpret as evidence of the gradual dispersal of mono-$\mathrm{[Fe/H]}$ populations, presumably due to radial migration and radial heating. The variation in shape of the broken exponential profile changes differently depending on the population $\mathrm{[Fe/H]}$, with low $\mathrm{[Fe/H]}$ populations inner profiles flattening faster, whereas the high $\mathrm{[Fe/H]}$ outer profiles flatten faster. We interpret this effect as tentative evidence for $\mathrm{[Fe/H]}$ dependent radial migration arising from pre-existing $\mathrm{[Fe/H]}$ gradients in the star forming disk. We showed that our results qualitatively reproduce those of Hayden et al. (2015), finding a skewed MDF that varies as a function of R.

• •

Flaring We find that flaring seems to be present in almost all mono-age populations, at differing levels. We have shown that the inverse flaring scale length $R_{\mathrm{flare}}^{-1}$ increases with age, meaning that the youngest populations flare most strongly. This finding appears inconsistent with that above, under the assumption that flaring is the result of radial migration. However, these results may be reconciled by invoking a more active accretion history in the early life of the disk, which could have suppressed flaring (e.g. Minchev et al., 2014b).

• •

The surface-mass density at $R_{0}$ We have measured the surface mass-density at the solar radius for each mono-age, mono-$\mathrm{[Fe/H]}$ population, finding a total surface mass density of $\Sigma_{R_{0},\text{tot}}=20.0_{-2.9}^{+2.4}\mathrm{(stat.)}_{-2.4}^{+5.0}\mathrm{(syst.)}\ \mathrm{M_{\odot}\ pc^{-2}}$. Before allowing for systematics, this value is less than current estimates (e.g. Bovy et al., 2012a; Flynn et al., 2006; McKee et al., 2015), however, the systematic uncertainties are large, mainly due to a mismatch between the $\log{g}$ scales in APOGEE and the PARSEC models, and as such, we find our value to be consistent within the uncertainties. The relative contribution of high to low $\mathrm{[}\alpha\mathrm{/Fe]}$ populations, $f_{\Sigma}$, is $18\%\pm 5\%$, which is consistent with existing measurements (e.g. Bland-Hawthorn & Gerhard, 2016).

• •

The $h_{Z}$ distribution at $R_{0}$ The shape of the mass-weighted $h_{Z}$ distribution found by this study is in good agreement with that of Bovy et al. (2012a), calling into question the existence of a vertical structural discontinuity in the Milky Way disk. The reconciliation of this finding with the discontinuity in chemical space (e.g. the bimodality in $\mathrm{[}\alpha\mathrm{/Fe]}$ at fixed $\mathrm{[Fe/H]}$: Hayden et al., 2015; Nidever et al., 2014) may shed new light on our understanding of the formation of the Galactic disk.

• •

The surface-mass density profile of the Milky Way We have found the combined (from mono-age, mono-$\mathrm{[Fe/H]}$ populations at low and high $\mathrm{[}\alpha\mathrm{/Fe]}$) surface-mass density weighted profiles of the Milky Way disk as a function of $\mathrm{[}\alpha\mathrm{/Fe]}$, age and $\mathrm{[Fe/H]}$, and found that the total surface density is also described by a broken exponential. We find that our results fail to determine the sign of the inner exponential to high significance out to $\sim 10$ kpc, but detect a turnover to a declining exponential, at high significance, thereafter. We find evidence of a radial mean age and $\mathrm{[Fe/H]}$ gradient driven by the changing dominant population as a function of radius. A detailed comparison of these findings with numerical simulations is necessary for a proper interpretation. Our finding of a decline in stellar density may be consistent with that found in other studies (e.g. Reylé et al., 2009; Sale et al., 2010), albeit at shorter radii.

These findings are strongly constraining to future theoretical work. With the recent (Lindegren et al., 2016) and future releases of Gaia data, and the ongoing APOGEE-2 survey (Sobeck et al., 2014), which will include an updated APOKASC sample (Pinsonneault et al., 2014), access to improved positions, abundances and age estimates is within reach. We again stress here that the age uncertainties in this data set can be as large as $40\%$, and so until more precise ages are attained for similarly sized samples, our conclusions must be considered under the caveat that the mono-age populations at old age are likely mixed to some extent. It will be possible to investigate this issue better once better ages for a larger sample are released by APOGEE and APOKASC (Pinsonneault et al., 2014).

Future studies of simulations which accurately track chemical evolution, gas and stellar dynamics, and the feedback processes which are dominant in galaxies will no doubt lead to a deeper insight into the physical processes leading to the present day structure of the Milky Way. The understanding of discontinuity in chemical space, namely the bimodality in $\mathrm{[}\alpha\mathrm{/Fe]}$ at fixed $\mathrm{[Fe/H]}$, and how this can be reconciled with the apparent structural continuity which we find here poses an interesting challenge to models of the formation of the Milky Way disk. By performing a first mapping of the 3D distribution of stellar populations as a function of age, metallicity, and [a/Fe], we hope that this work provides the kind of data needed for a comparison with numerical simulations that is unencumbered by the complexities associated with corrections for the survey selection function.

Acknowledgements

It is a pleasure to thank Marie Martig for helpful comments and discussion, and also for providing the catalogue of ages for APOGEE DR12 in digital format. We also thank Ivan Minchev, Misha Haywood and Paola Di Matteo for insightful discussion and comments which improved the clarity of the manuscript. The analysis and plots in the paper used iPython, and packages in the SciPy ecosystem (Jones et al., 2001; Hunter, 2007; Perez & Granger, 2007; van der Walt et al., 2011). JB received support from the Natural Sciences and Engineering Research Council of Canada, an Alfred P. Sloan Fellowship, and from the Simons Foundation. SM has been supported by the Premium Postdoctoral Research Program of the Hungarian Academy of Sciences, and by the Hungarian NKFI Grants K-119517 of the Hungarian National Research, Development and Innovation Office.

Funding for SDSS-III was provided by the Alfred P. Sloan Foundation, the Participating Institutions, the National Science Foundation, and the U.S. Department of Energy Office of Science. The SDSS-III web site is http://www.sdss3.org/.

Appendix A Density Fits

For completeness, we briefly discuss the quality of the fits performed with the method outlined in section 3. Figures 16 and 17 show the distance modulus distribution of the APOGEE data in each of our mono-age and mono-$\mathrm{[Fe/H]}$ bins (grey histograms) and the resulting distance modulus distribution when the best fit density model for each bin is run through the calculated effective selection function (which is the space in which models are fit in our procedure). The red line represents a a single-exponential fit to the radial and vertical spatial distribution and the black lines give the best fit broken-exponential density model (upon which we base our results). We show the single-exponential fit in order to demonstrate that, in most cases, this does not provide a good fit to the data, and that when a single exponential is a better fit, the broken-exponential density fit matches it.

Regarding Figure 16, which shows the low $\mathrm{[}\alpha\mathrm{/Fe]}$ sub-populations, it is clear that the black curve (broken exponential) represents a far better model for the data than the red curve (single exponential), in all mono-age, mono-$\mathrm{[Fe/H]}$ bins. While the black curve is not perfect in all cases, the peak of the distribution tends to lie at the correct $\mu$, whereas the red curve finds a peak at higher $\mu$ in most cases (due to the higher than necessary density at low Galactocentric radius in this model).

Figure 17 demonstrates the fits for the high $\mathrm{[}\alpha\mathrm{/Fe]}$ sub-populations. The greyed out panels reflect those with less than 30 stars, which we deem too noisy to render reliable fits. In many of the remaining panels, the red curve is similar or identical to the black, due to the fact that many of the high $\mathrm{[}\alpha\mathrm{/Fe]}$ populations are better described by single exponentials, and the broken exponential generally recovers this result. In most of the cases where the curves differ greatly, the red curve recovers the peak of the distribution better than the black - suggesting that breaks which were fit in the radial range we consider are artificial, and due to the noisy data in this regime. We discuss the broken exponential fits in the main text in order to make proper comparison with the low $\mathrm{[}\alpha\mathrm{/Fe]}$ sample, although it seems plausible that the single exponential model provides a better explanation of the data.

Appendix B The effect of uncertainties on trends with age

In order to demonstrate and characterise the effect that the age errors have on our interpretation of the trends between structural parameters and age, we created a mock data set from a set of randomly sampled density distributions. We created a mock data set with an input trend of $h_{Z}$ with age which increased monotonically with age from 0.2 to 1.2 kpc. Ages were assigned to each $h_{Z}$ population, sampling uniformly in bins of width 2 Gyr, to which we then added a random gaussian error of $40\%$, replicating the shifting of stars with different $h_{Z}$ into each age bin. In each bin, we sample a single exponential (a broken exponential with $R_{\mathrm{peak}}=0$) with scale length 8 kpc. This higher scale length is required to make the test computationally efficient, to produce realistic numbers of stars when selecting stars in APOGEE fields and does not impact on the results of this test. We assume no error on the stellar positions, in order to isolate the effect of the age errors.

While this test is a somewhat simplistic representation of the underlying processes, it serves as a good example of the effect of the age uncertainties which are expected in the present data. One example of its simplicity is the assignment of a single $h_{Z}$ to relatively wide bins in age. It seems logical to assume that if there is an age-$h_{Z}$ relation, then the change in $h_{Z}$ should be somewhat continuous with age. Our test assigns the same $h_{Z}$ to stars at bin edges (which should have $h_{Z}$ close to that of the bin-edge stars in the neighbouring bin). This may artificially increase the amount of blurring of the age-$h_{Z}$ trend. We also simplify the test by assuming that the only changing parameter is the scale height. Realistic structural parameters would change the relative number of stars within each bin observed by APOGEE (and considered in our test), and may change the level of contamination between bins. However, we assume that our simple approximation represents a ‘worst case’ scenario, where the mixing between bins is maximal.

We restrict the mock data to the APOGEE fields, simplifying the selection to a distance cut (assuming the selection fraction is 1 out to a distance which corresponds to $M_{\mathrm{H}}=-1.5$, assuming no extinction). We apply the method described in Section 3 to our mock data, fitting a broken exponential profile, and using the best fit solution to initiate an MCMC sampling of the posterior probability distribution. As in the main body of the paper, the reported parameter values reflect the median and $\sigma$ of one dimensional projections of the MCMC chain. The resulting age-$h_{Z}$ relation is shown in Figure 18. A clear trend is recovered between age and $h_{Z}$. The trend is still recovered at high age, regardless of the high level of mixing between bins, which increases the size of the error bars. The higher-scale height components are recovered by the analysis, but results are scattered around the input values, with large error bars. This serves to show that even in the face of large age uncertainties causing mixing between the adopted bins, our method is still able to recover the underlying trends of parameters with age.