Exploring the Stellar Age Distribution of the Milky Way Bulge Using APOGEE

Sten Hasselquist; Gail Zasowski; Diane K. Feuillet; Mathias Schultheis; David M. Nataf; Borja Anguiano; Rachael L. Beaton; Timothy C. Beers; Roger E. Cohen; Katia Cunha; José G. Fernández-Trincado; D. A. García-Hernández; Doug Geisler; Jon A. Holtzman; Jennifer Johnson; Richard R. Lane; Steven R. Majewski; Christian Moni Bidin; Christian Nitschelm; Alexandre Roman-Lopes; Ricardo Schiavon; Verne V. Smith; Jennifer Sobeck

doi:10.3847/1538-4357/abaeee

1. Introduction

The vast majority of disk galaxies in the local universe harbor an overdensity of light in their centers, commonly referred to as a "bulge." These bulges appear with a variety of structures and stellar populations, presumably resulting from different evolutionary processes: spheroidal "classical" bulges dominate in the most massive disk galaxies, while pseudobulges (mostly bars) are more common in Milky Way (MW)-mass galaxies (e.g., Fisher & Drory 2011). When these bulges—barred or not—form and how they evolve over time are still open questions; likely, these complex systems grew through some combination of accretion of stars that now reside in the center of the MW (e.g., Tumlinson 2010) and in situ star formation in a disk that later buckled into the structures observed today (see review in Athanassoula 2005).

The MW itself has a central, asymmetric, boxy overdensity of light first measured in integrated infrared (IR) photometry (e.g., Weiland et al. 1994; Dwek et al. 1995). Additional photometric and kinematic studies soon revealed that this boxy, peanut-shaped structure is likely a result of seeing a bar structure edge-on (e.g., Hammersley et al. 1994; Athanassoula 2005; McWilliam & Zoccali 2010; Nataf et al. 2010; Wegg & Gerhard 2013) and that the inner MW harbors a barred mass distribution with a semimajor axis ≃5 kpc (e.g., Wegg et al. 2015; Bovy et al. 2019) and major–minor axis ratio of 0.4 (Bovy et al. 2019). This barred central structure, which we will refer to as the "bulge," contains ∼50% of the MW's stellar mass (Licquia & Newman 2015), with ∼60% of that bulge mass residing in the bar structure and 40% in the inner disk (e.g., Portail et al. 2017). In order to understand the full picture of how our MW Galaxy formed and evolved, we must understand the star formation and chemical-enrichment histories of this critical region that contains the majority of stellar mass (e.g., Rich 2013; McWilliam 2016; Barbuy et al. 2018).

The stars that reside in the bulge span ∼2.5 dex in metallicity (nearly the full range observed across the MW); the metallicity distribution function (MDF) changes dramatically as a function of position in the inner Galaxy, resulting in an average negative metallicity gradient with height above the MW plane (e.g., Zoccali et al. 2008, 2017, 2018; Johnson et al. 2013; Rojas-Arriagada et al. 2014; Fragkoudi et al. 2018; García Pérez et al. 2018). In addition, the kinematics of the stellar populations are correlated with their chemistry, with the lower-metallicity stars on more spheroidal orbits and the higher-metallicity stars on more "bar-like" orbits (e.g., Hill et al. 2011; Ness et al. 2013, 2016b; Zasowski et al. 2016; Barbuy et al. 2018). Despite this spatially variant, broad MDF, the detailed elemental abundances of the bulge stars appear to be not only relatively homogeneous throughout the bulge, but also nearly identical to the chemical abundance pattern of thick-disk stars at the solar radius (e.g., Rojas-Arriagada et al. 2017; Haywood et al. 2018; Zasowski et al. 2019, but see e.g., Johnson et al. 2014). An exception to this similarity arises from the nonnegligible fraction of inner Galaxy stars originating from dissolved globular clusters (e.g., Fernández-Trincado et al. 2017; Schiavon et al. 2017).

These chemodynamical patterns suggest that the bulge formed many of its stars early, in a rapid star-formation event (e.g., McWilliam & Rich 1994; Fulbright et al. 2007; Johnson et al. 2012), and probably in a disk that later buckled into a boxy bar, giving rise to the metal-rich stars on radial bar-like orbits (e.g., Athanassoula 2016). However, major uncertainties remain regarding the timing of this buckling, the detailed structure of the disk before buckling (e.g., Fragkoudi et al. 2017), and the extent to which star formation has proceeded since this buckling event.

To answer these questions, we need ages for large numbers of bulge stars at all Galactocentric radii (R_cy) and distances from the plane ( $| {Z}_{\mathrm{GC}}|$ ), including in the midplane itself. Most age studies of the bulge to date have shown that all, or nearly all, of the stars appear to be old (>9–10 Gyr; e.g., Zoccali et al. 2003; Clarkson et al. 2011; Barbuy et al. 2018; Renzini et al. 2018). However, other groups found evidence for significant fractions of bulge stars with ages <8 Gyr (e.g., van Loon et al. 2003; Bensby et al. 2013, 2017; Catchpole et al. 2016), which are expected from some simulations, especially if only stars in the plane are considered (e.g., Ness et al. 2014). Such a large fraction of intermediate-age stars is consistent with some color magnitude diagram (CMD)-based studies (e.g., Holtzman et al. 1993; Haywood et al. 2016), but inconsistent with others (e.g., Clarkson et al. 2009; Gennaro et al. 2015; Surot et al. 2019). We refer the reader to Nataf (2016) and Barbuy et al. (2018) for recent reviews, and we note that these statements only apply to the larger bulge population beyond the innermost few hundred parsecs, where there is known to be recent and ongoing star formation (e.g., Morris & Serabyn 1996; Longmore et al. 2013).

In addition to the ambiguous presence of relatively young stars in the bulge, the spatial variations of the mean stellar age and the stellar age distribution are highly uncertain. The majority of bulge age studies have been limited to pencil-beam fields, typically not in the high-extinction midplane. Recent work from Bovy et al. (2019) argues that the bar appears distinctly older (and more metal-poor) than the inner disk, with a mean age of ≃8 Gyr (see also Sit and Ness 2020). In contrast, Wegg et al. (2019) find evidence for a metal-rich bar population younger than the local disk. We discuss both of these findings in the context of our analysis in Section 4.

Large-scale spectroscopic surveys continue to provide an ever-growing amount of chemical, kinematic, and, most recently, age information for stars across the MW, allowing us to explore numerous formation scenarios in different parts of the Galaxy. The Apache Point Observatory Galactic Evolution Experiment (APOGEE; Majewski et al. 2017) is a survey particularly well equipped to study the inner MW's evolution. Because APOGEE operates in the near-IR, the survey can observe stars in the MW center, even through the thick dust in the midplane of the Galaxy. In addition, the abundance of carbon and nitrogen features in its spectra allows for precise determination of carbon and nitrogen abundances, which in turn can be correlated with asteroseismology-derived masses and mapped to stellar ages.

Numerous studies have explored the MW disk in this way, for example, Masseron & Gilmore (2015), Martig et al. (2016b), Ness et al. (2016a), Mackereth et al. (2019), and Hasselquist et al. (2019b). Works that map ages onto APOGEE stars have usually relied on the exquisite APOKASC (Pinsonneault et al. 2014) and APOKASC-2 (Pinsonneault et al. 2018) stellar masses and ages as a training sample. However, the APOGEE stars observed in the bulge are typically much cooler and more luminous than those APOKASC stars with precise masses, meaning that computing ages for bulge stars is impossible without significant extrapolation.

In this paper, we use a new training set to compute ages for ∼46,000 stars in the MW Galaxy, and we analyze the results for ∼6000 stars in the bulge, nearly all of which are beyond the parameter space of the APOKASC-2 sample. We describe our data and bulge sample selection in Section 2, the spectral age information in Section 3.1, our training set in Section 3.2, and our implementation of The Cannon in Section 3.3. Our results—including mean age maps, the presence of young stars, and on-/off-bar differences—are presented in Section 4 and discussed in Section 5.

2. Observations and Sample Selection

Observations were taken as part of the APOGEE (Majewski et al. 2017), part of the fourth iteration of the Sloan Digital Sky Survey (SDSS-IV; Blanton et al. 2017). The APOGEE instruments are high-resolution, near-infrared spectrographs (Wilson et al. 2019) observing from both the Northern Hemisphere at Apache Point Observatory (APO) using the SDSS 2.5 m telescope (Gunn et al. 2006) and the Southern Hemisphere at Las Campanas Observatory (LCO) using the 2.5 m du Pont telescope (Bowen & Vaughan 1973). As of 2019 December, the dual APOGEE instruments have observed some 500,000 stars across the MW, targeting these stars with selections described in Zasowski et al. (2013, 2017), with updates to the targeting plan described in F. A. Santana et al. (2020, in preparation) and R. L. Beaton et al. (2020, in preparation).

Spectra are reduced as described in Nidever et al. (2015) and analyzed using the APOGEE Stellar Parameters and Chemical Abundance Pipeline (ASPCAP; García Pérez et al. 2016). A detailed analysis of the accuracy and precision of the stellar parameters and abundances can be found in Holtzman et al. (2018) and Jönsson et al. (2018). Our analysis uses results from the 16th Data Release (DR16) of the SDSS collaboration (Ahumada et al. 2020), which is the first data release containing data from the Southern instrument. Further explanations and assessments of this data release, including quantification of potential offsets between the Northern and Southern spectrographs, can be found in Jönsson et al. (2020).

This work focuses on the APOGEE stars that reside in the bulge, which we define using Galactic cylindrical coordinates, R_cy = $\sqrt{{X}_{\mathrm{GC}}^{2}+{Y}_{\mathrm{GC}}^{2}}$ and $| {Z}_{\mathrm{GC}}|$ . We transform the APOGEE stars into these coordinates based on their position in the sky and distance, after adopting a solar Galactocentric position of X_GC = −8.125 kpc (Gravity Collaboration et al. 2018) and Z_GC = 20.8 pc, as was done in Bovy et al. (2019). We adopt the stellar distances derived using astroNN²⁷ (Leung & Bovy 2019). These distances have uncertainties ∼15% at the location of the bulge, but both Queiroz et al. (2020) and Bovy et al. (2019) noticed that the DR16 astroNN distances are slightly underestimated at locations >5 kpc from the Sun. Therefore, we apply the same distance correction as derived and implemented by Bovy et al. (2019).

We show the resultant spatial distribution of our stellar sample in both the X_GC–Y_GC plane and R_cy–Z_GC plane in the left two columns of Figure 1. We define the bulge to be all stars with R_cy < 3.5 kpc and $| {Z}_{\mathrm{GC}}| \lt$ 1.5 kpc. Only stars that meet the following criteria are plotted in Figure 1:

1.
Signal-to-noise ratio (S/N) > 70 per pixel
2.
[Fe/H] > −0.5²⁸
3.
No STAR_BAD bit of ASPCAPBAD flag set
4.
$\mathrm{log}(g)$ < 3.3.

The third panel of Figure 1 shows the $\mathrm{log}(g)$ distribution of the stars that fall inside this bulge cut for each distance sample (blue histogram). The APOKASC sample, which has been used in the past as a training set for data-driven age-determination methods, is shown in red. As discussed further in Section 3, previous works that have employed The Cannon or similar techniques to derive ages of APOGEE stars use training sets, such as APOKASC (Pinsonneault et al. 2014), that have few stars with $\mathrm{log}(g)$ < 2.0 and almost no stars with $\mathrm{log}(g)$ < 1.5. Therefore, to derive ages for these "luminous giants" ( $\mathrm{log}(g)$ < 2.0) that make up the bulge APOGEE sample, we must use a different age-training sample.

3. Stellar Ages

The goal of this work is to derive ages for the luminous giants ( $\mathrm{log}(g)$ < 2.0) that primarily comprise the APOGEE bulge sample. We use The Cannon (Ness et al. 2015) to derive ages, but because of the vastly different $\mathrm{log}(g)$ distributions highlighted in Figure 1, we require a new training set to derive ages using this tool. In this section, we confirm there is age information in the APOGEE spectra for the luminous giants (Section 3.1), describe a new training set we will use to derive labels with The Cannon (Section 3.2), and discuss the application of The Cannon to our data set (Section 3.3).

3.1. Age Information in APOGEE Spectra

The age information in the APOGEE spectra of red giant stars primarily comes from carbon and nitrogen molecular features. This is because the birth [C/N] abundance of a star is affected by first dredge-up as it ascends the red giant branch (RGB). This dredge-up operates in such a way that the resultant [C/N] abundance ratio after the star has ascended the red giant branch depends on the mass of the progenitor, such that more massive stars have lower [C/N] abundances (e.g., Gratton et al. 2000; Martell et al. 2008; Salaris et al. 2015). Stellar models can then be invoked to obtain an age for an RGB star of a given mass and metallicity. There are now several works in the literature where this dependence has been exploited to interpret [C/N] abundance variations across the Galaxy as age variations (Masseron & Gilmore 2015; Hasselquist et al. 2019a), and even works that map ages directly onto stars observed by APOGEE (e.g., Martig et al. 2016b; Ness et al. 2016a; Mackereth et al. 2019).

These studies that derive ages for APOGEE stars all rely on the APOKASC and APOKASC-2 sample (Pinsonneault et al. 2014, 2018) as a training set. These stars have precise masses (≲10%) derived using asteroseismology. Ages are then inferred using stellar evolutionary models. However, this APOKASC sample is limited in the parameter space covered; specifically, it is lacking in stars with $\mathrm{log}(g)$ < 2.0 and nearly completely devoid of stars with $\mathrm{log}(g)$ < 1.5. Therefore, it is not an ideal training set for deriving ages for the APOGEE bulge stars, the majority of which have $\mathrm{log}(g)$ < 2.0 (as shown in the right panel of Figure 1), as it often results in ages derived from model extrapolation. This is summarized in the top row of Figure 2, where we show that, while the age coverage of the APOKASC-2 sample is quite good in the [C/N]–[Fe/H] plane for stars on the lower giant branch (2.6 < $\mathrm{log}(g)$ < 3.3), there are far fewer stars with $\mathrm{log}(g)$ < 2.0. The stars that are there are old, on average, and do not span the full range of [C/N]–[Fe/H] space covered by the bulge sample, which highlights the need for a different training set.

**Figure 2.** The [C/N]–[Fe/H] plane colored by mean $\mathrm{log}(\mathrm{age})$ . The top row shows this plane for the APOKASC-2 age-training sample, and the bottom row shows this plane for the "Feuillet" ages, described more in Section 3.2. The left column compares training sets for stars below the red clump, and the right column compares training sets for stars above the red clump, which is where the majority of the APOGEE bulge sample studied here is found.
Download figure:
Standard image High-resolution image

**Figure 2.** The [C/N]–[Fe/H] plane colored by mean $\mathrm{log}(\mathrm{age})$ . The top row shows this plane for the APOKASC-2 age-training sample, and the bottom row shows this plane for the "Feuillet" ages, described more in Section 3.2. The left column compares training sets for stars below the red clump, and the right column compares training sets for stars above the red clump, which is where the majority of the APOGEE bulge sample studied here is found.
Download figure:
Standard image High-resolution image

3.2. Luminous Star Training Set

To derive bulge ages using The Cannon, we must find a large sample of high-luminosity ( $\mathrm{log}(g)$ < 2.0) stars that have ages. Fortunately, several studies in the literature (e.g., Feuillet et al. 2016, 2018; Anders et al. 2018; Queiroz et al. 2018) have shown that it is possible to use Bayesian isochrone matching to derive precise (∼0.2 dex uncertainty in $\mathrm{log}(\mathrm{age})$ ) ages for red giant stars if the distance is known to better than 10%. With Gaia DR2 (Brown et al. 2018), there are now several thousand low-gravity stars in the APOGEE sample with distances more accurate than 10%. These stars are distributed across the MW disk over the range 6 kpc < R_cy < 11 kpc, resulting in a training set that spans a range of stellar populations. We refer to the ages derived in the methods described below as "Feuillet ages."

3.2.1. Deriving Ages for the Training Set

Briefly, ages for stars in the APOGEE luminous giant training set were derived using a simple Bayesian isochrone-matching method described by Jørgensen & Lindegren (2005), which produces an age probability distribution function (PDF). To derive the PDF for each star, we compute a likelihood function by comparing the measured parameters from each star with scaled-solar PARSEC (PAdova and TRieste Stellar Evolution Code) isochrones (Bressan et al. 2012) and an assumed prior star formation history (SFH), MDF, and initial mass function (IMF). The age of each star is assigned to be the age-weighted mean of the final PDF, and the age uncertainty is the formal dispersion in the PDF.

The measured parameters used for the likelihood function are effective temperature, surface gravity, α-adjusted metallicity, and absolute K magnitude. The effective temperature and surface gravity are taken as the calibrated values provided in APOGEE DR16. The α-adjusted metallicity is calculated based on the calibrated metallicity ([M/H]) and α-element abundance ([α/M]) from APOGEE DR16 using the formula of Salaris et al. (1993). The absolute K magnitude is calculated using the 2MASS K magnitude (Skrutskie et al. 2006), the distance from Bailer-Jones et al. (2018), and the K-band extinction provided by APOGEE, AK_TARG (or AK_WISE if unavailable).

The prior on the SFH is flat in $\mathrm{log}(\mathrm{age})$ as all calculations are done in $\mathrm{log}(\mathrm{age})$ , and the PARSEC isochrones used are gridded in $\mathrm{log}(\mathrm{age})$ with step sizes of 0.05. The MDF is assumed to be flat across the observational metallicity uncertainty, which is small compared to the typical spread of the disk MDF. We use the Chabrier lognormal IMF (Chabrier 2001) provided with the PARSEC isochrones. When analyzing a large sample of stars, a selection-function term is usually also included to account for the imposed limitations on surface gravity, color, and other parameters. However, the stars in the training sample in the current work are far from the edges of these selection cuts and are minimally affected. A full description of the age-determination method is described in Feuillet et al. (2016, 2018). Note we do not perform the hierarchical modeling step employed by Feuillet et al. (2016) for the present work.

3.2.2. Validation

The bottom row of Figure 2 shows how well the [C/N]–[Fe/H] plane is covered by the training set of Feuillet ages for both the lower giant branch (2.6 < $\mathrm{log}(g)$ < 3.3, left panel) that overlaps the APOKASC sample and the luminous giants ( $\mathrm{log}(g)$ < 2.0, right panel). The bottom left panel qualitatively appears similar to the upper left panel, indicating that ages are comparably mapped to the [C/N]–[Fe/H] space for both the APOKASC-2 seismic ages and the Feuillet ages for the lower giant branch stars. The right column emphasizes the lack of coverage using the APOKASC-2 ages as compared to the Feuillet ages, and also shows that the Feuillet ages map in a similar way to the [C/N]–[Fe/H] plane for both the lower and upper giant branches, confirming that there is age information in the C and N abundances for the luminous giants. Because the Feuillet ages rely on precise parallaxes of stars generally found in the MW disk, these two samples roughly sample the same volume of the Galaxy (inside ∼3 kpc from the Sun).

Throughout this work, we only analyze stars with [Fe/H] > −0.5. This is motivated by both theoretical studies of metallicity-dependent extra mixing along the giant branch (e.g., Carbon et al. 1982; Charbonnel & Lagarde 2010) and recent empirical measurements of this extra mixing in the APOGEE data by Shetrone et al. (2019). Shetrone et al. (2019) show that stars with $\mathrm{log}(g)$ < 2.0 and [Fe/H] < −0.5 can have extra mixing further affect the [C/N] abundance ratio by ∼0.15 dex at [Fe/H] = −0.7 and 0.58 dex at [Fe/H] = −1.4, making the [C/N] abundance ratio less sensitive to age. While we do not use the [C/N] abundance ratios explicitly to derive ages in this work, we know much of the age information comes from C and N features, likely making our ages susceptible to the same reduction in sensitivity.

An additional complication of deriving ages using carbon and nitrogen spectral features or [C/N] abundances is the difficulty of taking into account potential birth abundance variations in [C/N] abundances across the Galaxy. The APOKASC sample is confined to one region of the disk, and the Feuillet ages come from nearby stars (within ∼3 kpc from the Sun) by design, where the birth [C/N] abundance does not seem to vary by much (see, e.g., Martig et al. 2016a). We discuss how a varying birth [C/N] in the bulge might affect our results in Section 5.4.1.

3.3. The Cannon with Feuillet Ages

Having verified that age information is contained in the APOGEE spectra for the luminous giants and that our training set spans a parameter space similar to the APOGEE bulge sample, we can use The Cannon (Ness et al. 2015) to derive ages for these bulge stars. The Cannon is a data-driven technique for deriving stellar labels, where a model describing the flux at each pixel is created from a training set with well-known stellar labels. This model is then applied to a "test" set of spectra where the derived labels are returned. In this work, we fit a quadratic model to the spectra in our training set. While the [C/N] abundance could be used to map ages directly onto APOGEE stars via a multiparameter fit, as was done in Martig et al. (2016a), we opt to use the entire spectral range, as Ness et al. (2016a) showed that mass/age information is also encoded in the ¹²C/¹³C ratio.

3.3.1. Method

We build an input training sample from all stars with Feuillet ages (described in Section 3.2) using the following cuts:

1.
S/N > 100 per pixel
2.
[M/H] > −0.5
3.
No STAR_BAD bit of ASPCAPFLAG flag set
4.
No STARFLAG bits set
5.
0.5 < $\mathrm{log}(g)$ < 2.0
6.
Gaia DR2 parallax uncertainty <10%.

These cuts result in a training sample of 3711 stars that span the parameter space shown in Figure 3. We use T_eff, $\mathrm{log}(g)$ , [M/H], [Mg/Fe], and $\mathrm{log}(\mathrm{age})$ as the input labels. While the parameter space is reasonably well covered, we note here that the ages will be extrapolated for stars with [M/H] > +0.4 as the training set contains very few stars this metal rich. Therefore, ages for stars with [M/H] > +0.4 should be used with caution. We run The Cannon using code obtained from Anna Ho,²⁹ which also renormalizes the APOGEE spectra. This additional normalization step changes the APOGEE-normalized spectra very little, but it does ensure that the spectra obtained from the Northern and Southern instruments are normalized in the same way.

We first validate the output of the model by training The Cannon using 90% of this training sample, then deriving labels for the remaining 10%. We do this 10 times and analyze how the input labels compare to the output labels from this cross-validation test. The results are shown in Figure 4. Pearson correlation coefficients and standard deviations of the differences between the input and output labels are shown in the upper left of each panel.

**Figure 4.** Density plots showing the results of the 90–10 cross-validation test for each label. Pearson correlation coefficients and standard deviations are indicated in the upper left of each panel.
Download figure:
Standard image High-resolution image

We find that we are able to reproduce T_eff, $\mathrm{log}(g)$ , [M/H], and [Mg/Fe] to high precision, but the resultant ages are less precise, with a scatter around the 1–1 line of ∼0.3 dex, implying an uncertainty of ∼0.22 dex, which is slightly higher than other age studies using similar methods. However, these are the first such ages where the training set sufficiently covers the bright end of the RGB, so that ages for our bulge stars are not extrapolations of the method. We further explore the age accuracy and precision in Section 3.3.2 and the Appendix.

Using this training set, we run The Cannon on ∼46,000 luminous red giant stars in the APOGEE sample. These stars are selected to cover the same range of label space as the training set, specifically as follows:

1.
S/N > 70 per pixel
2.
0.5 < $\mathrm{log}(g)$ < 2.0
3.
−0.5 < [Fe/H] < +0.5
4.
No BAD bit of ASPCAPBAD flag set.

When deriving labels for the test sample, we additionally remove any stars from the training sample that had reduced χ² > 2 in the cross-validation step, and we retrain the model. Ages for all ∼46,000 stars along with the DR16 [Fe/H] measurements and uncertainties are provided in Table 1. The full table can be found in machine-readable format in the online journal.

Table 1. Ages Produced in This Work

APOGEE ID	APOGEE Field	Telescope	$\mathrm{log}(\mathrm{age})$	σ_log(age)	[Fe/H]	σ_[Fe/H]
2M00000002+7417074	120+12	apo25m	9.36	0.30	−0.17	0.01
2M00000317+5821383	116-04	apo25m	9.68	0.28	−0.28	0.01
2M00000546+6152107	116+00	apo25m	9.02	0.24	−0.27	0.01
⋯^a	⋯	⋯	⋯	⋯	⋯	⋯

Note.

^aThe full table can be found in machine-readable format on the online journal.

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as: Data Typeset image

3.3.2. Validation

The 90–10 test conducted in the previous section on the training sample suggests that we are reliably recovering age information from the APOGEE spectra. However, there are several external checks we can do to assess our age precision and accuracy further. Additional details can be found in the Appendix, but we summarize the results here.

First we check for potential age dependence on telescope and instrument setup. Using a sample of 62 stars that were observed in both hemispheres for which we were able to derive ages, we find that, while most labels are identical, the age labels are offset such that the stars observed from LCO are 0.08 dex younger in $\mathrm{log}(\mathrm{age})$ than the APO stars. Therefore, we apply a 0.08 dex offset to all LCO stars. We apply the offset in this direction because ∼75% of our training set is composed of stars observed from APO. This offset has already been applied to the ages reported in Table 1, and we note here that this offset does not affect our conclusions significantly.

The results of the cross-validation test (Figure 4) suggest a precision in $\mathrm{log}(\mathrm{age})$ of ∼0.22 dex. After exploring the potential dependence of this precision on S/N, [M/H], [Mg/Fe], and $\mathrm{log}(g)$ , we find that the precision only depends on $\mathrm{log}(g)$ , with the lower- $\mathrm{log}(g)$ stars being less precise. In Appendix A.2 we describe how we derive a quadratic fit to $\mathrm{log}(\mathrm{age})$ random uncertainties as a function of $\mathrm{log}(g)$ . These uncertainty values are reported in Table 1 and range from 0.3 dex precision for the lowest $\mathrm{log}(g)$ values ( $\mathrm{log}(g)$ ≃ 0.5–0.8) to 0.2 dex precision for the highest $\mathrm{log}(g)$ values ( $\mathrm{log}(g)$ = 2.0).

For a first check of the accuracy of our age results, we cross-match our age catalog to the open cluster catalog of Cantat-Gaudin et al. (2018). We find three clusters for which we have derived ages for more than three members: NGC 6791, NGC 6819, and NGC 2204. We find that the ages we derive using the median of the ages of cluster members agree to ∼0.1 dex or better for the youngest two clusters, but we find a median age that is 0.2 dex higher than the age reported in the Kharchenko et al. (2013) catalog for NGC 6791. These results are explained in more detail in Appendix A.3.

We can also check the accuracy of our results by comparing stars that have ages both from The Cannon and APOKASC, limited to 2.0 < $\mathrm{log}(g)$ < 1.5. We find that the ages agree reasonably well where the APOKASC mass uncertainties are low, but when the APOKASC mass uncertainties become larger than ∼10%, the APOKASC ages are systematically larger than the ages we derive using The Cannon. This is discussed in more detail in Appendix A.4, but given that the offset is only for stars with large APOKASC uncertainties and on the upper end of the $\mathrm{log}(g)$ range considered, we do not apply any offset to our data.

For another check on the accuracy of our ages, we reproduce age-abundance maps of the MW that have been previously studied in works such as Ness et al. (2016a) and Martig et al. (2016b). These maps are shown and described in further detail in Appendix A.5. Although we study more luminous giants with ages derived from a different training set than all previous work, we find qualitative trends very similar to that first found and described by Ness et al. (2016a) and Martig et al. (2016b). We also subselect our sample to match the spatial distribution of the CoRoGEE sample studied by Anders et al. (2017) and measure near-identical metallicity gradients for stars younger than 2 Gyr and stars older than 10 Gyr.

Finally, in Appendix A.6, we show the ages we are able to get for the most metal-rich stars in the Large Magellanic Cloud (LMC) and Sagittarius Dwarf (Sgr) galaxies. We find that the most metal-rich LMC stars have a median age of ∼0.8 Gyr, consistent with the expected ages of these stars from star formation history studies. Similarly, we find the Sgr stars to have age ∼6 Gyr, again consistent with what is expected from the literature.

4. Results

We now present our age results, first showing mean age maps for much of the entire APOGEE sample, then focusing on the bulge sample defined in Section 2.

4.1. Mean Age Maps

First we introduce some qualitative spatial age trends of the bulge. In Figure 5 we show maps of the X_GC–Y_GC plane of the MW for different $| {Z}_{\mathrm{GC}}|$ bins colored by number density of stars (Σ, top row), [Fe/H] (second row), $\mathrm{log}(\mathrm{age})$ (third row), and η (fourth row), which is the number of standard deviations each bin i is away from the mean $| {Z}_{\mathrm{GC}}|$ of that entire $| {Z}_{\mathrm{GC}}|$ range that divides the columns. This serves as a metric to assess potential spatial biases induced in the binning:

$\begin{eqnarray}&&\eta =\displaystyle \frac{\langle {Z}_{\mathrm{GC}}{\rangle }_{i}-\langle {Z}_{\mathrm{GC}}{\rangle }_{\mathrm{Tot}}}{{\sigma }_{| {Z}_{\mathrm{GC}}| }}.\end{eqnarray} \tag{ 1 }$

The top row of Figure 5 shows the stellar density at each position in the Milky Way to demonstrate how APOGEE samples the bulge at these metallicities and $\mathrm{log}(g)$ values (also see Bovy et al. 2019 and Queiroz et al. 2020). As expected, the near side of the bulge (X_GC < 0 kpc) is covered at a much higher stellar density than the far side (X_GC > 0 kpc). The far side also has relatively incomplete X_GC and Y_GC coverage, especially for the two bins with $| {Z}_{\mathrm{GC}}|$ < 0.5 kpc. This is not surprising given the extinction in the midplane, where A_V frequently exceeds 25 mag (e.g., Schultheis et al. 1999). Therefore, for all future analysis discussed in this work, we restrict our sample to be on the near side of the bulge, or X_GC < 0 kpc.

The second row of Figure 5 shows how the mean [Fe/H] of our sample, which is limited to [Fe/H] > −0.5, changes with position in the Galaxy. As has been found before for stars at R_cy < 8 kpc, the mean [Fe/H] decreases with increasing $| {Z}_{\mathrm{GC}}|$ (e.g., Hayden et al. 2015). Similar to Leung & Bovy (2019), we find that the metallicity of the MW appears to peak around R_cy ∼ 4–5 kpc, and maybe even decreases in the innermost region. We also see that the region of the Galaxy with 0 kpc < Y_GC < 3 kpc and −3 kpc < X_GC < 0 kpc appears to be more metal-poor than other regions of the bulge. This low-metallicity feature was also seen by Leung & Bovy (2019) and A. B. A. Queiroz et al. (2020, in preparation). Bovy et al. (2019) interpreted this as a signature of the bar, which also was shown to be distinct in age and kinematics. However, Wegg et al. (2019) actually find the bar to be more metal-rich than the disk. This is potentially a result of incomplete $| {Z}_{\mathrm{GC}}|$ sampling, described more below.

The third row of Figure 5 shows how the mean $\mathrm{log}(\mathrm{age})$ changes with position in the Galaxy. Stars outside 3 kpc appear to exhibit a steeper vertical age gradient than stars inside 3 kpc, which appear to all be old. This is in qualitative agreement with studies finding the mean age of bulge stars with [Fe/H] > −0.5 to be old. This results in a radial age gradient where the ages of stars in the plane ( $| {Z}_{\mathrm{GC}}| \,\lt$ 0.25 kpc) go from $\mathrm{log}(\mathrm{age})$ ∼ 9.8–9.9 at R_cy < 3 kpc to $\mathrm{log}(\mathrm{age})$ ∼ 9.3 at R_cy = 5 kpc. Out of the plane, there appears to be no radial age gradient in 0 kpc < R_cy < 8 kpc. In the planar bin, $| {Z}_{\mathrm{GC}}|$ < 0.25 kpc, we also find signs of the off-bar side of the bulge (l < 0°, Y_GC < 0 kpc) being younger, on average, than stars on the on-bar side of the bulge (l > 0°, Y_GC < 0 kpc), as was found by Bovy et al. (2019). We further explore and quantify these differences in Section 4.4.

However, the interpretation of such gradients, as well as on-bar versus off-bar comparisons, is complicated by potential selection biases. We show in the bottom row of Figure 5 the same maps but colored by η, defined as the number of σ(Z_GC) values away from the mean of the entire $| {Z}_{\mathrm{GC}}|$ range (Equation (1)). The X_GC–Y_GC bins that have little color have stars with a mean $| {Z}_{\mathrm{GC}}|$ close to the mean, whereas bins that are dark red or dark blue have stars that are up to 2σ away from the mean. In the case for the $| {Z}_{\mathrm{GC}}|$ < 0.25 kpc bin, we show that the bulge is probed at a larger mean $| {Z}_{\mathrm{GC}}|$ than the disk region (R_cy > 4.0 kpc), and different sides of the bulge are probed at different mean $| {Z}_{\mathrm{GC}}|$ heights. Therefore, any intrinsic vertical age/metallicity gradients that exist in the inner Galaxy will potentially cause one to measure different mean ages/metallicities if the $| {Z}_{\mathrm{GC}}|$ sampling is not taken into account. We discuss how this affects spatial variations in the age distribution of the bulge in Section 4.4.

4.2. Age–Metallicity Relation of the Bulge

For the remainder of this section, we will focus on the bulge stars only, which we define as follows:

1.
R_cy < 3.5 kpc
2.
$| {Z}_{\mathrm{GC}}|$ < 1.5 kpc
3.
X_GC < 0 kpc (near side of MW only).

The reason that we restrict our sample to the near side of the bulge is the "patchy" nature of the APOGEE bulge coverage, as shown in Figure 5. Figure 6 shows the age–metallicity relation for these stars in the same $| {Z}_{\mathrm{GC}}|$ bins shown in Figure 5. In the upper panel, we calculate the linear metal enrichment over time, Δ_Z = dZ/dt, where Z here is the metallicity mass fraction, not Galactic height, Z_GC. To calculate this value, we find the mean metallicity of stars with 9.75 < $\mathrm{log}(\mathrm{age})$ < 9.95 in each Z_GC bin and assume the stars are enriched from metallicity Z = 0 at t = 13.7 Gyr to the Z observed at this age range, which is ∼7 Gyr ago. We choose this age to compare to other works that typically calculate this enrichment over this time period. This age limit is usually imposed because after 7 Gyr, the stars enrich very slowly over time (e.g., Bernard et al. 2018) or are not even found in other samples.

We find that the metallicity evolution of the bulge was quicker for the stars closer to the plane than for stars out of the plane. We also find that our stars with Z_GC < 0.5 kpc have Δ_Z values consistent with the model put forth by Haywood et al. (2016), but shallower than the fit for the same age range in Bernard et al. (2018).

We also overplot the Bensby et al. (2017) stars in the middle two bins, corresponding to the most likely spatial overlap given the latitude of the Bensby et al. (2017) stars. The ages of these stars are not inconsistent with what we find in our sample. We also see that the supersolar metallicity stars extend down to the lower ages that Bensby et al. observed. As described more in Section 4.3, we find that the mean age of the stars with +0.2 < [Fe/H] < +0.4 is ∼9.5–9.6 in $\mathrm{log}(\mathrm{age})$ , or 3–4 Gyr old. This is only true for the spatial bin closest to the plane.

We also find signs that the ages of the most metal-poor stars (−0.5 < [Fe/H] < −0.3) in our sample become slightly younger, on average, at larger distances from the plane. Bernard et al. (2018) find a reasonable spread in age at these metallicities, suggesting that finding stars at −0.5 < [Fe/H]< −0.3 that are as young as 7–8 Gyr old is not unusual. The fact that these stars are younger farther from the plane also fits with an overall slower metallicity evolution farther from the plane and represents the metallicity of the last stars formed at these $| {Z}_{\mathrm{GC}}|$ heights.

4.3. Young Stars in the Bulge

The age–metallicity relations suggest that the bulge is not uniformly old, and we find stars as young as ∼1–3 Gyr at [Fe/H]> 0.1 and $| {Z}_{\mathrm{GC}}| \,\lt$ 0.25 kpc. To assess the significance of these stars, given our age uncertainties, we divide our sample into monoabundance bins (e.g., Bovy et al. 2016), and we analyze the age distribution of each monoabundance bin. Specifically, we are interested in the median age of each monoabundance bin and whether each monoabundance bin can be described by a single age.

Figure 7 shows the median age (top row), standard deviation in age (second row), and skewness in age (third row) for each monoabundance bin for the four $| {Z}_{\mathrm{GC}}|$ bins. Stars are divided into 0.1 dex bins of [Fe/H] for −0.5 < [Fe/H] < +0.5, and 0.2 dex bins of [Mg/Fe] for −0.1 < [Mg/Fe] < +0.5. Points are plotted according to the median [Fe/H] and [Mg/Fe] of each bin, so they do not necessarily lie on a grid point. Points are only plotted if the bin contains more than 20 stars.

As shown in the top row of Figure 7, we find the bins containing stars with [Fe/H] > +0.2 and [Mg/Fe] < +0.1 are generally the youngest in median age, with the age of these populations slightly increasing from in the plane ( $\mathrm{log}(\mathrm{age})$ ∼ 9.5) to out of the plane ( $\mathrm{log}(\mathrm{age})$ ∼ 9.8). While stars outside these abundances are generally all older with $\mathrm{log}(\mathrm{age})$ > 9.9, there is potentially a slight age gradient, such that stars with [Fe/H] < −0.2 are ∼0.1 dex younger out of the plane than stars in the plane.

We also show the standard deviation of $\mathrm{log}(\mathrm{age})$ for each monoabundance population in the middle row of Figure 7. Given our age uncertainties of 0.22–0.30 dex, points that are colored black are consistent with being composed of stars of the same abundance and age. Points that are lighter in color likely have a real intrinsic age spread. We find nearly no age spread for monoabundance bins at large distances from the plane, but substantial age dispersion for monoabundance bins in the plane, particularly at [Fe/H] < 0.0 and supersolar [Fe/H] stars with [Mg/Fe] < +0.1. In the plane, the populations with [Fe/H] > +0.2 and [Mg/Fe] < +0.1, which are the populations with the youngest median age, actually exhibit a smaller age dispersion than the more metal-poor stars also with [Mg/Fe] < +0.1.

Finally, the third row of Figure 7 shows the same monoabundance bins, but now colored by the skewness of the age distribution. Dark blue points correspond to abundance bins that contain a skew toward younger ages, and dark red points correspond to abundance bins that contain a skew toward older ages. The two youngest bins in the plane actually exhibit a slight positive skew, indicating that these abundance bins still contain old stars. The stars at large distances from the plane with [Fe/H] < −0.1 exhibit a large negative skew, suggesting that these more metal-poor stars are not uniformly old and contain a smaller fraction of younger stars, as also found by Bernard et al. (2018).

In Figure 8 we summarize the $| {Z}_{\mathrm{GC}}|$ heights and [Fe/H], where excess young stars can be found in the bulge. We plot the fraction of stars younger than 8 and 5 Gyr, as was done in Bernard et al. (2018), but do this for a range of $| {Z}_{\mathrm{GC}}|$ heights. All $| {Z}_{\mathrm{GC}}|$ bins exhibit an increase in the fraction of younger stars with increasing [Fe/H], but this increase starts at lower [Fe/H] for stars closer to the plane. This is true for both 8 and 5 Gyr. The gray lines in each panel of Figure 8 show the expected fraction, given our age uncertainties, if the stars in each [Fe/H] bin were formed at a single age of 9, 10, and 11 Gyr, with a 0.1 dex spread in $\mathrm{log}(\mathrm{age})$ . Also, the left panel of Figure 8 shows that stars with $| {Z}_{\mathrm{GC}}|$ > 0.5 kpc and [Fe/H] < −0.3 have a nonnegligible fraction of stars with age <8 Gyr, potentially a consequence of the overall slower chemical evolution farther from the plane.

So while many of the stars, especially at solar [Fe/H] and below, appear to be consistent with being born from one event some 9–10 Gyr ago, we measure a significant fraction of younger stars for the more metal-rich stars, where "young" means age <5 Gyr. However, we note that the youngest stars in our full MW age sample at these metallicities are actually found outside the bulge (R_cy ∼ 5 kpc, see Appendix A.5). Therefore, while we see little evidence for very recent star formation in the bulge, we do find strong evidence that the bulge formed stars as recently as 2–5 Gyr ago, and that this star formation took place at +0.2 < [Fe/H] < +0.4 and in the plane.

To further analyze the spatial and kinematic properties of the younger stars we find in the bulge, we divide the bulge sample into 0.1 dex bins of [Fe/H], and we define an "old" and "young" sample for each bin, where "old" stars are stars greater than one standard deviation away from the mean age of each bin, and "young" stars are stars less than one standard deviation away from the mean of each bin. The median R_cy, $| {Z}_{\mathrm{GC}}|$ , and v_ϕ of these stars are plotted as a function of the median [Fe/H] in Figure 9.

The first row of Figure 9 shows that the young stars with +0.1 < [Fe/H] < +0.3 tend to be found at larger Galactocentric radii than the old stars at the same metallicity and at slightly lower $| {Z}_{\mathrm{GC}}|$ (second row of Figure 9). The young stars at these metallicities also tend to have larger rotational velocities than the older stars (third row of Figure 9), suggesting that the young stars are generally still in a disk as compared to the old stars at those metallicities. Again, we find that the younger metal-poor stars are typically found at higher $| {Z}_{\mathrm{GC}}|$ than the older stars.

To summarize, while we do indeed find that many stars in the bulge at the metallicities studied here are 9–10 Gyr old, there is a statistically significant number of metal-rich stars younger than 5 Gyr that tend to be at 2 < R_cy < 3 kpc and $| {Z}_{\mathrm{GC}}| \,\lt$ 0.25 kpc and exhibit kinematics consistent with rotating around the MW center. Therefore, while the bulk of the metal-rich bulge formed in one event some 9–10 Gyr ago, star formation continued in a disk at supersolar metallicities until ∼2 Gyr ago.

4.4. On and Off the Bar

Because APOGEE-2 has observed stars across the entire bulge, we are able to compare stars on the on-bar side of the bulge to those on the off-bar side. To select these stars, we follow the lead of Bovy et al. (2019) and put stars in the on-bar sample if they fit inside an ellipsoidal structure oriented at 25° from the Sun, with a half-width of 2 kpc. The off-bar sample consists of the stars that fall outside this region. This ellipse used to denote on- and off-bar stars is shown on the maps of Figure 5. To account for potential radial variation, we also limit both samples to be 2.0 < R_cy < 3.5 kpc and, as before, only consider stars at X_GC < 0 kpc.

The mean maps shown in Figure 5 suggest that the off-bar side is more metal-rich and slightly younger than the on-bar side, but only at $| {Z}_{\mathrm{GC}}|$ < 0.50 kpc. We measure these differences and quantify their significance in Figure 10. The top row shows the age distribution of the stars in the off-bar (red) and on-bar (blue) samples defined above. While the median ages of the two samples are nearly identical, the results of a Kolmogorov–Smirnov (K-S) test suggest a small but potentially significant difference in the age distributions. The differences become more pronounced when we limit the sample to $| {Z}_{\mathrm{GC}}|$ < 0.50 kpc (middle row of Figure 10). However, the right panel highlights that there is also a difference in the $| {Z}_{\mathrm{GC}}|$ distributions probed, leaving open the possibility of sampling biases.

**Figure 10.** Age distributions and $| {Z}_{\mathrm{GC}}|$ distributions shown separately for the on-bar (blue) and off-bar (red) samples in bins of $| {Z}_{\mathrm{GC}}|$ . Results of K-S tests are indicated in the upper left of each panel. Median $\mathrm{log}(\mathrm{age})$ values are included in the left column panels.
Download figure:
Standard image High-resolution image

We correct for this by subsampling the on-bar sample, which contains ∼10 times more stars than the off-bar sample, such that it matches the same $| {Z}_{\mathrm{GC}}|$ distribution as the on-bar sample. The resultant age distributions are shown in the bottom row of Figure 10. With this correction, we find that the K-S test cannot rule out the possibility that these two distributions are drawn from the same parent distribution. More data from the Southern Hemisphere telescope will help to further quantify the similarities (or differences) of these on- and off-bar age distributions, as well as more careful separation into on- and off-bar groups using orbital information (see A. B. A. Queiroz et al. 2020, in preparation).

The above results also suggest that the apparent old age of the bar shown in Figure 5 is largely or entirely a sampling effect, with the stars in the inner 0.1 kpc of the plane exhibiting younger ages, on average, than the rest of the stars. Additionally, because the mean age of the on-bar sample decreases when we subsample to match the $| {Z}_{\mathrm{GC}}|$ distribution of the off-bar sample, then there is likely at least a slight vertical age gradient. Quantifying this vertical age gradient is beyond the scope of this work, as careful evaluation of selection effects would need to be taken into account.

5. Discussion

5.1. Spatial Age Trends of the Bulge

We find that ∼50% of our bulge sample, composed of stars with −0.5 < [Fe/H] < +0.5, are older than 8 Gyr. As shown in Figure 8, nearly all stars with subsolar metallicity have ages consistent with being born from one star formation event some 9–10 Gyr ago. This agrees with many lines of evidence that point to an old bulge, including the chemical track of bulge stars in [α/Fe]–[Fe/H] space (e.g., Cunha & Smith 2006; Fulbright et al. 2007; Hill et al. 2011; Johnson et al. 2011; McWilliam 2016; Bensby et al. 2017; Bovy et al. 2019; Zasowski et al. 2019), where most or all of the subsolar-metallicity stars are shown to be enhanced in their α-elements. Recent work suggests that the chemistry of the bulge stars is identical to the chemistry of the thick-disk stars found in the Solar Neighborhood (e.g., Rojas-Arriagada et al. 2014; Bovy et al. 2019). These stars are found to be at ages similar to what we find for our bulge stars at subsolar metallicities, serving as further evidence for a coeval formation scenario.

We do find some spatial variations in the age distributions of the bulge. First, we find a slight vertical age gradient at 0 < $| {Z}_{\mathrm{GC}}|$ < 0.5 kpc, where the stars closer to the plane are ∼0.2 dex younger than stars farther from the plane. A much steeper vertical gradient is observed outside the bulge at R_cy ∼ 5 kpc, which is driven by the inclusion of the youngest stars (0.5–2 Gyr) in the plane. These youngest stars are largely not present in the bulge sample. We find that the bulge is enriched at a rate of dZ/dt ∼ 0.0034 Gyr⁻¹, which is very similar to the rate predicted by Haywood et al. (2018), but we find that this rate decreases with height above the plane.

The older mean age and spatial age variation we observe in this work is qualitatively similar to that recently found by Bovy et al. (2019). Their work likely uses stars similar to what we use here, but the ages are derived in a very different fashion, with different training sets. Bovy et al. (2019) describes the bar as a structure that is distinct in mean age. This is also a prediction of Ness et al. (2014), who find that the mean age of the off-bar stars in the plane is younger than the mean age of the on-bar stars in the plane. Our maps displayed in Figure 5 show results similar to Bovy et al. (2019), in that spatial bins corresponding to the on-bar side of the MW are older than for the off-bar. However, as described in Section 4.4, once we correct for different $| {Z}_{\mathrm{GC}}|$ sampling, we cannot conclude that the on-bar and off-bar populations are different in age at $| {Z}_{\mathrm{GC}}|$ < 0.25 kpc. Future data from sight lines at b = 0° of the bulge combined with more careful selections (see A. B. A. Queiroz et al. 2020, in preparation) will more conclusively answer this question.

In addition to a vertical age gradient in the bulge, Ness et al. (2014) find that the youngest stars should be in the plane at $| {Z}_{\mathrm{GC}}|$ < 0.14 kpc. We do find a nonnegligible fraction of young stars that share the following properties:

1.
+0.2 < [Fe/H] < +0.4
2.
2.0 kpc ≲ R_cy ≲ 3.0 kpc
3.
$| {Z}_{\mathrm{GC}}|$ < 0.5 kpc
4.
v_ϕ ∼ 150 km s⁻¹, or ∼50 km s⁻¹ larger than older stars at the same metallicity.

Therefore, our picture of the bulge from this work is that stars across all metallicities were formed some 8–10 Gyr ago and can be found all across the bulge. The more planar regions of the bulge were enriched more quickly than the off-plane regions of the bulge, with the outermost $| {Z}_{\mathrm{GC}}|$ bin exhibiting no stars with [Fe/H] ⪆ +0.2. However, superposed on this old stellar population is a younger, temporally extended population that appears to be in the plane and at supersolar metallicities. So, after the bulge/bar formed in the initial burst, star formation still occurred, at a lower rate, in a disk or ring, until stopping some 2–4 Gyr ago.

5.2. Young Star Reconciliation

The paradigm for many decades was that the bulge only contained old stars, where "old" refers to stars with age ≳8 Gyr. As eloquently summarized in Nataf (2016), multiple lines of evidence exist for an exclusively old bulge, as well as multiple lines of evidence for some young-to-intermediate-age, metal-rich bulge stars. Because we find our younger stars in the plane to be more metal-rich and rotating around the Galactic Center, then this means that studies that are biased against metal-rich stars, look at sight lines at b > 5°, or remove stars based on high proper motions may be removing these younger stars. Given these criteria, we try to reconcile the confirmed presence of younger stars in this work, and other works such as Bensby et al. (2013), with other studies that find no younger stars.

Kuijken & Rich (2002) reanalyzed the Baade's Window Hubble Space Telescope data of Holtzman et al. (1993), whose initial conclusions suggested the supersolar metallicity stars have a median age of ∼5 Gyr. Baade's Window is at b ∼ −4°, so there should be some younger, metal-rich stars in these fields. However, in their reanalysis, Kuijken & Rich (2002) found no young stars after cleaning this sample with new proper-motion measurements. While a detailed sight line-by-sight line proper-motion analysis is beyond the scope of this paper, the young stars we identify would have reasonably large proper motions in the l direction at the longitude of Baade's Window, suggesting a strict proper-motion cut could potentially remove these stars.

A similar argument could be made for why young stars are not found in the Zoccali et al. (2003) and Clarkson et al. (2011) samples, although the former sample appears to contain fewer supersolar-metallicity stars—likely a consequence of studying a higher latitude field, where the metal-rich stars are a weaker component of the bulge MDF (e.g., Figure 6, Zoccali et al. 2017; García Pérez et al. 2018). Future analysis, consisting of a field-by-field proper-motion comparison, will determine whether or not the proper-motion cleaning is a reason why these photometric studies are lacking in younger bulge stars. However, Haywood et al. (2016) suggest that some of these photometric studies (e.g., Valenti et al. 2013), which show a tight main-sequence turnoff, actually necessitate a decent fraction of young, metal-rich stars in the bulge (see also Barbuy et al. 2018).

If the works that show an exclusively old bulge can be explained by probing sight lines higher from the plane or proper-motion cuts inducing a bias, then the works that show younger stars better fit into the paradigm listed above. The Bensby et al. (2013) sample consists of stars that are nearly all at $| b|$ < 5°, so they are likely at a $| {Z}_{\mathrm{GC}}|$ height where we see young stars, but we do not know the exact distance to these stars. As already shown in Figure 6, the younger stars in both of our samples overlap. Younger metal-rich stars are also found by Bernard et al. (2018), who do use proper-motion cuts, but they study stars at 2 $\lt \,| b|$ < 4°. The presence of younger and intermediate-age stars at latitudes closer to the plane is also consistent with the findings of Catchpole et al. (2016), although they find that the long-period Mira variables (age ∼5 Gyr) exhibit a clumpy distribution and suggest that they are associated with the bar.

In summary, the likelihood of observing younger stars in the bulge depends on the metallicity and $| {Z}_{\mathrm{GC}}|$ probed. However, nowhere do we find stars as young as we do in the disk outside of the bulge region (see Appendix A.5). Therefore, in the areas where we find younger bulge stars, star formation was still shut off some ∼2+ Gyr ago, so any analyses that find statistically significant numbers of stars in the bulge with age <2 Gyr would still be difficult to reconcile with the present analysis, as well as most or all other studies in the literature.

5.3. The Milky Way in the Galactic Context

While we would agree that the bulge region of the MW contains numerous old stars, as most studies of external galaxies find, we show that it does depend on the sample location and what metallicities are being probed. Kruk et al. (2018) find that, after decomposing the inner regions of their barred galaxy sample into bar+disk, the disk component is often bluer or younger. This agrees qualitatively with what we see in the MW. Additionally, Fragkoudi et al. (2020) find that some of the galaxies studied in their barred galaxy sample from Auriga simulations exhibit "inner rings" of recent star formation, resulting in a nonnegligible fraction of stars with ages <5 Gyr.

There have also been many studies to try to understand how the barred regions of external galaxies differ from the disk regions in both age and metallicity gradients. Some works show that the gradients along bars are shallower than gradients off the bar (e.g., Fraser-McKelvie et al. 2019; Neumann et al. 2020), while other works do not show this (e.g., Sánchez-Blázquez et al. 2014). See Seidel et al. (2016) for a more complete description of these discrepancies. A detailed age-gradient study is beyond the scope of this work, but we do find some evidence suggesting a vertical age gradient in at least the off-bar region of the MW, and possibly in the on-bar region as well (see Figure 5 and Section 4.4). Perhaps this can resolve the tension on whether or not one observes radial age-gradient variations for on-bar and off-bar populations in external galaxies, as the result will depend on distance from the plane probed.

Finally, the fact that we observe a slower metallicity evolution with increasing height suggests that, while a bar may be efficient at washing out radial gradients, some $| {Z}_{\mathrm{GC}}|$ gradient still exists in the SFH that has remained measurable at present times.

5.4. Primordial Carbon, Nitrogen, and Helium Abundance Variations

Although we do not explicitly map [C/N] to ages in this work, The Cannon mainly cues off of C and N spectral features to derive the ages of stars (see, e.g., Ness et al. 2016a). The ages we derive then implicitly rely on the fact that C and N are dredged up and mixed in similar ways for the bulge stars as they are for stars near the solar circle, which are the stars that make up our training set.

Therefore, three potential caveats of this study could affect our age determination:

1.
The birth abundance of C and N for the bulge stars differs from the birth abundance of C and N for the solar circle stars, ultimately resulting in a different post-first-dredge-up [C/N] abundance.
2.
The birth abundance helium affects the way C and N are dredged up.
3.
The birth abundance helium itself varies enough to affect the ages of the stars in a drastic way, which we are not taking into account.

It is most important to consider the effect these caveats would have on our identification of younger stars in the bulge.

5.4.1. Birth Abundance C and N Variation

It is possible that the mass-dependent dredge-up that allows for the mapping of [C/N] to mass (and then age) may be affected by large differences in the birth C and N abundances, or even more specifically, the birth [C/N] abundance ratio. Martig et al. (2016a) used a sample of APOGEE subgiants, which span no more than 2–3 kpc from the solar circle, to show that the birth [C/N] abundance ratio did not vary wildly with Galactic position. Unfortunately, no large sample of subgiants with precise C and N abundances exists for the bulge.

Despite this, one argument we can make for minimal birth [C/N] abundance variation relative to our training-set [C/N] abundances is that the bulge chemical abundances appear to be very similar to the chemical abundances of the high-α or thick-disk stars that are found within the solar circle. These stars appear to share very similar chemical abundance tracks in the α-elements, light odd-Z elements, and Fe-peak elements like Ni and Mn (e.g., Rojas-Arriagada et al. 2017; Zasowski et al. 2019; A. B. A. Queiroz et al. 2020, in preparation). It would be difficult to form stars that have the same α-element abundances but vastly different C abundances, or stars that have similar metallicity-dependent odd-Z element abundance patterns but vastly different N abundance patterns. Therefore, because the bulge sample is made up of stars with α-element abundances and metallicities that are represented in our training set, we have no reason to think that the birth abundance [C/N] varies in a way that would cause erroneous ages.

5.4.2. Helium Affecting Dredge-up

There is some work in the literature that suggests that helium can affect the dredge-up in giant stars. Karakas & Lattanzio (2014) find that, during the third dredge-up, an asymptotic giant branch star that is enhanced in helium will dredge up less carbon. However, we are most concerned about after the first dredge-up, as the majority of our stars should be RGB stars given the log(g) range studied. Salaris et al. (2015) found little to no effect of helium abundance on the [C/N] abundance after the first dredge-up, but they only considered small changes in helium (ΔY ∼ 0.02). There is some evidence that the bulge contains helium-enhanced stellar populations (e.g., Buell 2013; Nataf et al. 2013), so this work should be repeated with ΔY values of ∼0.06.

To investigate the potential effects of helium abundance on [C/N] abundance after the first dredge-up, we analyze stellar models from Tayar et al. (2017). The results are shown in Figure 11.

In the left panel of Figure 11, we show how the [C/N] abundance ratio changes as a 1.3 M_⊙, [Fe/H] = 0.2 star ascends the red giant branch for two different helium birth abundances. The star that is enhanced in helium ends up with a lower [C/N] abundance ratio after the first dredge-up. This means The Cannon might assign a younger age than is warranted. However, as shown in the text in the left panel of Figure 11, the star itself has a younger age because it is helium-enhanced.

To investigate which effect wins out, we show the ages of stars across a range of masses at [Fe/H] = 0.2 in the middle panel of Figure 11 for solar helium (blue) and enhanced helium (red) as a function of [C/N] after the first dredge-up. Both helium abundances follow a similar log(age)–[C/N] relation, with a slight offset between the two such that at fixed [C/N] the helium-enhanced stars are actually slightly younger. We fit a quadratic function to the solar-helium stars and derive a log(age)–[C/N] relation. We use this relation to infer an age of the helium-enhanced stars based on their [C/N] after first dredge-up. The right panel of Figure 11 shows the difference between the true age of these helium-enhanced stars and the ages inferred based on [C/N]. We find that the inferred ages are ∼0.06 dex older than the true ages. Therefore, our age method would actually put helium-enhanced stars slightly older than they actually are.

5.4.3. Helium Affecting Stellar Age

Stars enhanced in helium live shorter lifetimes. Bazan & Mathews (1990) suggest that the lifetime of a star can change by ∼0.3 dex in $\mathrm{log}(\mathrm{age})$ for helium Y = 0.3 for lower-mass stars. If the [C/N] abundance ratio after first dredge-up is not significantly impacted by birth helium abundance, then our age method would be pushing younger, helium-enhanced stars to older ages. This would mean that our sample would contain even more younger stars than what we currently find.

After all this, even if the helium conspires to give us younger ages in one way or the other, then instead of young stars, we have simply found "helium-peculiar" stars, and they tend to be found at +0.2 < [Fe/H] < +0.4, low $| {Z}_{\mathrm{GC}}|$ , R_cy = 2–3 kpc, and so on.

6. Summary

We summarize our conclusions below:

1.
We have derived ages for 47,000 luminous giant stars across much of the MW disk and bulge.
2.
The stars at −0.5 < [Fe/H] < +0.5 that reside in the bulge are primarily old, with about half of the stars in our sample being older than 8 Gyr.
3.
We find that the stars in the plane are enriched at a rate of dZ/dt ∼ 0.0034 Gyr⁻¹, which is exactly what was found by Haywood et al. (2016). However, stars out of the plane are enriched at a much slower rate (dZ/dt ∼ 0.0013 Gyr⁻¹, a factor of three slower), suggesting an inside-out bulge formation scenario where the stellar populations have not yet been fully mixed.
4.
We find a nonnegligible fraction of younger stars (2–5 Gyr old) that primarily have +0.2 < [Fe/H] < +0.4, $| {Z}_{\mathrm{GC}}| \,\lt$ 0.25 kpc, 2.0 kpc ≲ R_cy ≲ 3.5 kpc, and kinematics more consistent with rotation. This suggests an extended star-formation history of the bulge that took place in a disk after the initial burst.
5.
This work suggests that some of the literature disputes regarding an old versus young bulge can be settled by understanding that the younger stars are only found at $| {Z}_{\mathrm{GC}}|$ ≲ 0.5 kpc and at [Fe/H] ≳ 0.2.
6.
When correcting for $| {Z}_{\mathrm{GC}}|$ sampling effects, we do not find a measurable age difference between stars inside the bar and stars outside the bar, but future data will improve this measurement, as well as measurements of the vertical and radial gradients.

We thank the anonymous referee for their comments that improved this manuscript. We thank Melissa Ness for the extremely helpful guidance on how best to run The Cannon. We thank Jamie Tayar for the help in providing and explaining her stellar models.

S.H. is supported by an NSF Astronomy and Astrophysics Postdoctoral Fellowship under award AST-1801940. D.K.F. was supported by the grant 2016-03412 from the Swedish Research Council. D.A.G.H. acknowledges support from the State Research Agency (AEI) of the Spanish Ministry of Science, Innovation and Universities (MCIU) and the European Regional Development Fund (FEDER) under grant AYA2017-88254-P. A.R.L. acknowledges substantial financial support for this research from FONDECYT Regular 1170476 and QUIMAL project 130001. T.C.B. acknowledges partial support from Grant PHY 14-30152 (Physics Frontier Center/JINA-CEE), awarded by the U.S. National Science Foundation. D.G. gratefully acknowledges support from the Chilean Centro de Excelencia en Astrofísica y Tecnologías Afines (CATA) BASAL grant AFB-170002. D.G. also acknowledges financial support from the Dirección de Investigación y Desarrollo de la Universidad de La Serena through the Programa de Incentivo a la Investigación de Académicos (PIA-DIDULS). D.M.N. acknowledges support from NASA under award No. 80NSSC19K0589 and support from the Allan C. and Dorothy H. Davis Fellowship.

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

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

This research made use of Astropy,³⁰ a community-developed core Python package for Astronomy (Astropy Collaboration et al. 2013; Price-Whelan et al. 2018), SciPy (Jones et al. 2001), NumPy (Oliphant 2006), and Matplotlib (Hunter 2007).

Appendix: Age Validation

A.1. North versus South

While tests of the two APOGEE spectrographs suggest that their performance is nearly identical (Wilson et al. 2019), there are small variations in the line-spread function (LSF) across the detectors in both instruments, as well as variations in the LSF between the two instruments. Therefore, we might expect some differences in our results for Northern spectra than for Southern spectra, especially since the training set is ∼75% Northern spectra. Because several stars we derive ages for were observed from both the Northern and Southern instrument setups, we are able to quantify potential systematic uncertainties in age derivation based on whether the star is observed from the Northern or Southern Hemisphere.

These comparisons are shown in Figure A1. Only those stars with S/N > 70 from both hemispheres are plotted. The T_eff, $\mathrm{log}(g)$ , and [M/H] labels output by The Cannon are nearly identical. There may be a slight bias in derived [Mg/Fe], such that the APO [Mg/Fe] values are 0.02 dex lower than the LCO values. The ASPCAP [Mg/Fe] values for the same set of stars differ by 0.01 dex such that the APO values again are ∼0.01 dex lower than the LCO values. For the analysis in this paper, we use the ASPCAP [Mg/Fe] values, but note that we only use these values to bin stars into monoabundance bins, and adopting either set of [Mg/Fe] values has no significant effect on our results.

There is also a slight bias in $\mathrm{log}(\mathrm{age})$ , such that the APO ages are 0.08 dex older than the LCO ages, on average. In our analysis, we apply a 0.08 dex offset to the LCO ages to bring this overlap sample into better agreement. This offset is included in the ages provided in Table 1. This offset does not significantly influence any of our results.

A.2. Deriving Uncertainties

The 90–10 cross-validation test done in Section 3.3.1 suggests an age precision of ∼0.22 dex (standard deviation of the differences, 0.31, divided by $\sqrt{2}$ ). However, we find that this precision is not uniform across parameter space. In particular, the precision varies strongly with $\mathrm{log}(g)$ . We approximate this variation by calculating the standard deviation of the differences between the input age label and the label obtained during the 90–10 cross-validation step. We calculate this for 20 star bins in a "moving boxcar" fashion, and we fit a function to these standard deviations as a function of $\mathrm{log}(g)$ . The distribution of differences between two samples is related to the uncertainty of a single sample by $\sqrt{2}$ , so these standard deviations are divided by $\sqrt{2}$ . We fit the following function to these results, and this is the prescription for the uncertainties reported in Table 1:

$\begin{eqnarray*}&&{\sigma }_{\mathrm{log}(\mathrm{age})}=0.42\mbox{--}0.19\mathrm{log}(g)+0.035{\left(\mathrm{log}(g)\right)}^{2}.\end{eqnarray*}$

A.3. Comparison to Open Clusters

We match our age sample to the Gaia cluster catalog provided by Cantat-Gaudin et al. (2018). We find three clusters for which we have derived ages for more than three stars: NGC 6791 (nine stars), NGC 6819 (six stars), and NGC 2204 (six stars). We only select stars that have Cantat-Gaudin et al. (2018) membership probabilities of >0.5. We compare the median age we derive for the stars in these clusters to the ages provided in the Kharchenko et al. (2013) cluster catalog. The log(age) values we adopt from this catalog are 9.65 for NGC 6791, 9.21 ± 0.024 for NGC 6819, and 9.29 ± 0.018 for NGC 2204. The log(age) values we derive in this work, including the dispersion in log(age) for these clusters, are 9.86 ± 0.23 for NGC 6791, 9.11 ± 0.12 for NGC 6819, and 9.13 ± 0.07 for NGC 2204.

We find that we slightly underestimate the ages for the youngest two clusters by ∼0.1 dex, but we overestimate the age of the oldest cluster, NGC 6791, by ∼0.2 dex. However, some literature work suggests that NGC 6791 could be as old as log(age) ∼ 9.9 if helium enhancement is accounted for (e.g., Brogaard et al. 2012), which would agree much better with the age we would find for NGC 6791 using our nine members.

A.4. Comparison to APOKASC

We compare our ages to those derived from asteroseismic masses in the APOKASC-2 sample (Pinsonneault et al. 2018). This overlap sample is relatively small (412 stars) and only contains stars with $\mathrm{log}(g)$ between 1.5 and 2.0, but it serves as one of the only external checks we have for accuracy. Figure A2 shows that the ages agree well within the scatter for stars with APOKASC-2 mass uncertainty <10%, and there is no apparent systematic offset. However, above an APOKASC-2 mass uncertainty of 10%, there is an offset such that the APOKASC-2 ages are on average older than The Cannon ages. A detailed understanding of this apparent offset with uncertainty is beyond the scope of this work, but other studies that have not removed APOKASC-2 stars with large mass uncertainties have found that the ages they reproduce are generally younger than the input APOKASC-2 ages at older ages (e.g., Martig et al. 2016a; Mackereth et al. 2019), and these works typically apply their own correction factors.

**Figure A2.** Differences between APOKASC-2 $\mathrm{log}(\mathrm{age})$ and The Cannon $\mathrm{log}(\mathrm{age})$ ages derived in this work, plotted as a function of fractional APOKASC-2 mass uncertainty.
Download figure:
Standard image High-resolution image

We do not apply a correction to bring our results into better agreement with the APOKASC ages because it is the APOKASC stars with large mass uncertainties that are discrepant, and the vast majority of stars studied in this paper have $\mathrm{log}(g)$ < 1.5, so this potential offset is only motivated for a small fraction of our sample.

A.5. Reproducing MW Disk Results

To show further support for the accuracy of the ages derived in this study, we reproduce age maps of the MW that have been produced in other works. Our sample is different in that the ages are derived using a different training set than these other works, and we study the more luminous giants, whereas these other works typically study red clump giants or giants farther down the giant branch. These maps are shown in Figure A3. As mentioned in Section 4, the youngest stars in our full sample are not found in the bulge, but at R_cy > 3 kpc and $| {Z}_{\mathrm{GC}}| \,\lt$ 0.25 kpc. As has been found in other studies, the youngest stars in each radial bin in the planar region of the Galaxy (bottom row of Figure A3) are found at lower [Fe/H] from the inner to outer Galaxy (e.g., Martig et al. 2016b; Ness et al. 2016a; Mackereth et al. 2019). We also select stars in the same Galactic region as the CoRoGEE sample and measure a metallicity gradient of stars with age <2 Gyr and a gradient of stars with age >10 Gyr. We find these gradients to be −0.068 ± 0.001 dex kpc⁻¹ and −0.030 ± 0.003 dex kpc⁻¹, respectively, which are both in excellent agreement with the same gradients measured by Anders et al. (2017).

**Figure A3.** [Mg/Fe]–[Fe/H] maps of the MW colored by log(age). Columns are separated by bins of R_cy, and rows are separated by bins of Z_GC.
Download figure:
Standard image High-resolution image

We also see that stars increase in age from near the plane to above the plane, but the stars above the plane themselves exhibit a radial age gradient, qualitatively similar to what was found by Martig et al. (2016b).

A.6. LMC and Sgr Ages

APOGEE has also observed stars in the LMC and Sagittarius dwarf galaxy (Sgr), allowing us to further validate our ages. Because we are limited to deriving ages for stars with [Fe/H] > −0.5, we can only look at the ages for the most metal-rich stars in these two galaxies. The most metal-rich stars of these galaxies have [Fe/H] ∼ 0.0, so the following analysis is limited to stars with −0.5 < [Fe/H] < 0.0.

In Figure A4 we compare the age distributions of the two dwarf galaxies (LMC in red and Sgr in blue) to the age distribution of a low-latitude MW disk sample (green) and a high-latitude MW disk sample (orange). LMC stars are selected from the Nidever et al. (2020) sample, and Sgr stars are selected from the Hayes et al. (2020) sample. The MW stars are split in latitude to select a "thin disk" sample (∣b∣ < 4°) and a "thick disk" sample (∣b∣ > 20°).

We find that the LMC stars are young, with a median age of ∼0.8 Gyr. This is consistent with the expected age of these most metal-rich LMC stars based on the star formation history of Harris & Zaritsky (2009), as well as the the chemical evolution model invoked in Nidever et al. (2020) to explain the observed [Mg/Fe]–[Fe/H] abundance pattern. The Sgr stars are older, with a median age of ∼6 Gyr. This is consistent with N-body simulations of Sgr, which generally find that Sgr fell into the MW some 5–6 Gyr ago (e.g., Law & Majewski 2010). Moreover, the SFH from Siegel et al. (2007) suggests that the bulk of Sgr stars at these metallicities should be 4–8 Gyr old. We also find that, as expected, the thin-disk MW stars are younger than the thick-disk MW stars (∼2 Gyr old as compared to ∼6 Gyr old).