The Flattening Metallicity Gradient in the Milky Way's Thin Disk

Ages are calculated as in Vickers & Smith (2018), using a method similar to Jørgensen & Lindegren (2005). A grid of isochrones is constructed in age, metallicity, and α abundance. The grid is spaced in increments of 0.25 Gyr from 1 Gyr to 5 Gyr, and 0.5 Gyr thereafter to a maximum of 14.5 Gyr; 0.1 dex in metallicity from −2.4 to +0.5 dex; and 0.2 dex in α abundances from −0.2 to +0.8 (for [Fe/H] > 0 dex, [α/Fe] is limited to the range of −0.2 to +0.2).

For this work we have opted to use the Dartmouth library of isochrones (Dotter et al. 2007, 2008) because it has implemented α abundances and also has Gaia color filters. One downside of the Dartmouth family of isochrones is that it does not extend to very young stellar ages, with the youngest isochrone being 1 Gyr.

Each star is compared to each isochrone, where a χ² distance is calculated between the star and each line segment on the isochrone with respect to the observational errors. This χ² value is converted into a likelihood, which is then weighted by the initial mass function of Kroupa (2001) corresponding to each line segment of the isochrone. The integral of these likelihood values over the full mass range of the isochrone gives the likelihood of the star belonging to that specific isochrone. We then marginalize over all [Fe/H] and [α/Fe] values to find the posterior probability of a star being a given age (Equations (1)–(9)) in Vickers & Smith 2018, with [α/Fe] values accounted for in addition to [Fe/H]).

When comparing magnitudes to isochrones, we use Gaia photometry that has been extinction corrected with the mwdust ⁴ library detailed in Bovy et al. (2016). This library determines extinction values at any given coordinate using the 3D reddening maps of Drimmel et al. (2003), Green et al. (2015), and Marshall et al. (2006).

We use distances and errors from the catalog of Bailer-Jones et al. (2021). The asymmetric errors in the distance catalog are accounted for by sampling a two-piece Gaussian, that is, a broken Gaussian with different σ values, but the same probability mass, on each side of the mean. This preserves the probability mass on both sides of the catalog median but suffers from discontinuity at the center. An alternative option that we did not use is the "Fechner" distribution, which again uses different σ values to describe the positive and negative halves of the distribution, but stitches the two Gaussian halves together at the peak. The "Fechner" distribution preserves continuity at the mean but does not preserve the balance of probability mass to the positive and negative directions.

To increase the speed of the calculation, isochrone points are only considered when they are within 10 times the error of the stellar parameter in question.

Errors on our age calculations are derived following the procedure outlined in Jørgensen & Lindegren (2005). We first calculate the likelihood of a star being each age in our grid of possible ages (0.25, 0.5 Gyr, etc.), we then linearly interpolate the likelihoods between each of these points to create a likelihood surface. Next, we find the maximal extent in age such that all likelihoods outside of that extent are less than 0.61 times the maximum likelihood (there may be some likelihoods within that extent that are below 0.61 times the maximum likelihood, as in the case where a star is well described by two different age estimates). This cut of 0.61 has been found to well describe the 68% confidence intervals by Jørgensen & Lindegren (2005).

For those unfamiliar with stellar ages, there are a few things to keep in mind about isochrone stellar ages. The first is that areas of the parameter space where isochrones change dramatically with age (for example, the turnoff or subgiant branch) generally give more precise ages than areas that are similar for all ages (such as the cool main sequence or the giant branch); this is why some authors carefully select turnoff and subgiant samples when working with stellar ages. The second is that observed atmospheric parameter errors will naturally cause stars to be scattered outside of the extents of the isochrone grid, which leads to pileups of stars in the edge age bins. The third would be that younger isochrones change faster than older isochrones, which means younger stars generally have smaller absolute age errors than older stars (which is why many authors quote age errors as a percentage of age, rather than a Gyr value, or work with logarithmic ages). Some of these quirks will be referenced throughout the paper.

To investigate the accuracy and possible biases in our ages, we compare our ages to those of open clusters. We first select all of our data that have been assigned cluster membership in the catalog of Cantat-Gaudin et al. (2018) and then compare the ages of these cluster members to the literature ages of the clusters given in the catalog of Dias et al. (2021). There are 17 clusters in our comparison that have more than one member in our data and are present in both these catalogs.

In Figure 1 we show the distribution of clusters in metallicity space and the difference in age between our calculation and the literature values. While there appears to be an offset between our ages and the literature values given for the clusters (ours being about 0.66 Gyr older), there appears to be no significant trend between age offset and metallicity.

Zoom In Zoom Out Reset image size

It is relatively well known that the Dartmouth family of isochrones "runs hot" or, alternatively "runs old." We find that our Dartmouth isochrone ages are about 1.8 Gyr older than the ages derived by Sanders & Das (2018), who use the PARSEC isochrones (see Section 2.1.3), for example.

Some aspects of the test in Section 2.1.1 prompt closer investigation. For example, in Figure 1, our cluster with the most members and also one of the largest age offsets is NGC 2682 (M67). We see an offset of 3.24 Gyr between the cluster's literature value (3.76 Gyr) and the ages we calculate for the constituent stars (7.0 Gyr). However, the literature value is perhaps not consistent with more recent studies of the cluster, such as that of Barnes et al. (2016), which places the cluster's age closer to 4.2 Gyr and thus our offset closer to 2.8 Gyr.

We further explore the offset in this particular cluster in Figure 2. We plot our 90 member stars, color coded by their age determinations, along with three isochrones that all have [Fe/H] = −0.1 and [α/Fe] = 0.0, the closest isochrone in our grid to the member stars' values. We see that the literature-age isochrone does not align with the majority of the stars in the lower turnoff and main-sequence region. An intermediate-age isochrone overlays the data well, except for the upper turnoff region. An old-age isochrone misses the turnoff entirely.

Zoom In Zoom Out Reset image size

A pattern is visible in the individual stars, with log(g) appearing correlated with age, although comparison with the overplotted isochrones makes this correlation appear not unreasonable. Note that this figure does not show absolute magnitude, which is another constraint in the age-fitting procedure.

It is common in stellar age studies to encounter an age–metallicity degeneracy where younger ages can be misinterpreted as lower metallicities (as both parameters tend to move stars blueward or toward higher temperatures). This degeneracy could possibly explain the offset we see in age (with our ages being systematically older than the literature values) and [Fe/H] (with our metallicities being systemically lower than the literature values).

Examining our [Fe/H] data more closely, we see that there is a trend in metallicity, in that our data tend to grow more metal poor, and younger, as we travel up the isochrone from the main sequence to the turnoff. This trend is actually the opposite of the expected age-[Fe/H] degeneracy pattern.

These contradictory patterns—that our stars are too old and too metal poor for the cluster as a whole, and yet as we increase the offset in metallicity within the cluster, the age difference grows smaller—prevent us from making a simple correction to these data. The topic deserves more attention than we give it here, but unfortunately, it seems to be beyond the scope of the current work.

While our data has ages for this cluster that do not match the literature value, we see no obvious reason to distrust our age calculations in Figure 2. Rather, we reiterate to the readers that stellar ages rely on many different complex systems working in tandem, and findings relying on them should be considered in that light.

The determination of stellar ages is a challenging problem with many complexities (Soderblom 2010). For the interested readers, we briefly compare our age determinations with those from two other major catalogs: the LAMOST catalog of Xiang et al. (2017), and the composite age catalog of Sanders & Das (2018). Owing to the various quality cuts, procedures, and data sources used in these three different works, we have a sample of about 200,000 stars with ages in all three catalogs (rather than the ∼million or more in any given one of them).

In Figure 3 we plot the age versus age and age difference versus age plots for the data in these age catalogs, and we see general agreement between all three. Our catalog agrees best with the ages derived by Sanders & Das (2018), with a fairly constant offset of about 1.8 Gyr (where our ages are older). There is a more significant and nonuniform difference between our derived ages and those of Xiang et al. (2017), which may in part arise from those authors' unique approach to calibrating the stellar temperature values (to the scale of Huang et al. 2015; see their paper for more details), as well as their usage of different isochrones in their age calculation.

Zoom In Zoom Out Reset image size

This comparison is particularly interesting as these three works used three different sets of isochrones. As mentioned, we use the Dartmouth family of isochrones. Xiang et al. (2017) utilize the Yonsei–Yale isochrones (Yi et al. 2001; Kim et al. 2002; Yi et al. 2003; Demarque et al. 2004), and Sanders & Das (2018) make use of the PARSEC isochrones (Bressan et al.2012).

Various differences exist between these families of stellar models. α abundances, for example, are available for Dartmouth and Yonsei–Yale isochrones but are not yet readily available for Padova isochrones. Young stellar ages (younger than 1 Gyr) are available for Yonsei–Yale and PARSEC isochrones, but not for Dartmouth models (although the MIST isochrones from the lead author of the Dartmouth isochrones should be adding α-element abundances soon). Dartmouth models are also generally found to "run hot," leading ages derived with Dartmouth isochrones to usually be slightly older than those derived with other models.

While a more detailed comparison is beyond the scope of this work, it is encouraging to see that these three catalogs produce similar age determinations, and it is interesting to see the magnitude of difference that can be expected from using different procedures and models in the determination of stellar age.

Moving beyond comparison with other age catalogs derived using isochrones, we may also compare our derived ages with those calculated using astroseismology. Astroseismic ages (ages calculated based on the pulsations of stars) are considered to be some of the most accurate age estimates possible with current data. The downside of astroseismic ages is that high-quality, time-series photometry is required, which limits the sample size considerably compared to other methods.

Luckily, LAMOST has observed a large number of objects in the Kepler field, which is where the majority of astroseismic age-work has been done. In Figure 4 we compare our derived ages with the astroseismic ages available in the second APOGEE-Kepler astroseismic catalog, APOKASC (Pinsonneault et al. 2018). The APOKASC catalog is mostly giant stars.

Zoom In Zoom Out Reset image size

We see reasonable qualitative agreement between the two catalogs at intermediate and older ages (as determined in this work), and a fair amount of disagreement among the giant stars that we calculate to be young. The disagreement is to be expected because isochrone age determinations are notoriously less reliable in areas where isochrones of different ages overlap, such as the giant branch or main sequence, and the APOKASC catalog contains mainly red giants.

While there is little to be done regarding the young-age disagreement with our current data set, we note that only 15% of our data are giants (having surface gravities less than 3). We also rerun the main analysis of this work with a subsample of data that avoids the giant branch and main sequence, in Section 4.1.4, and find our results with both the full sample and high-quality subsample to be similar.

We use a right-handed coordinate system with the Sun located along the x-axis at (X, Y, Z) = (−8.27, 0, 0) kpc. The solar system rotates with a velocity of (V^R , V^ϕ , V^Z) = (0,−236, 0) km s⁻¹ (alternatively, (U, V, W) = (0, 236, 0) km s⁻¹), and it has a motion relative to the local standard of rest of (U, V, W) = (13.0, 12.24, 7.24) km s⁻¹. These values are adopted from Schönrich & Aumer (2017).

We assign Cartesian coordinates for our stars using the on-sky coordinates and distances from the catalog of Bailer-Jones et al. (2021). Cartesian velocities are calculated following the procedure of Johnson & Soderblom (1987). From these Cartesian velocities, we calculate the spherical and cylindrical velocities in the usual way.

To account for measurement errors, we create 100 mock observations per star with the relevant parameters pulled from their error distributions and accounting for the correlations between the measurements (as reported in the Gaia astrometric catalog). The final phase-space coordinates and errors are the mean and standard deviation of these 100 Monte-Carlo error-expanded mock observations.

We calculate the orbits based on these final phase-space coordinates with the Python package galpy (described in Bovy 2015. ⁵ ) Within this package, we use the provided potential MWPotential2014, which is a superposition of an exponentially cutoff power-law bulge (McMillan 2011), a Miyamoto–Nagai disk (Miyamoto & Nagai 1975), and an NFW (Navarro–Frenk–White) halo (Navarro et al. 1997). The potential is modified to reflect our coordinate system values for solar position and velocity, noted above. We will also make use of the potentials of Irrgang et al. (2013) and McMillan (2017), which are available in the same code library and are similarly modified for our solar values.

For this paper, we will be measuring the iron abundance [Fe/H] as a function of the guiding radius, R_g, defined as the radius of a circular orbit with the same angular momentum as the observation within the same Galactic potential. This is slightly more precise than approximating R_g as the average of the orbital apocenter and pericenter radii and accounts for the effect whereby stars are generally closer to pericenter than apocenter along their orbits. This R_g-based gradient is less affected by blurring (epicyclic motions) than an R_present-based gradient, and its use can help reveal the intrinsic underlying gradient.

In our thin disk selection, we included a stringent cut on Z_max to select stars whose orbits are confined to within 250 pc from the plane. As we have mentioned, the mechanics of churning may favor stars on colder orbits, specifically stars with low vertical velocity dispersions (Vera-Ciro et al. 2014) or low eccentricities, therefore we suspect that Z_max should be strongly correlated with observational evidence of that. That is, we expect the stars on the vertically coldest orbits to show the strongest evidence for the [Fe/H] gradient flattening relative to their gradients at birth. This flattening over time could even cause us to observe older stars with shallower metallicity gradients than younger stars, despite presumably having been born with steeper radial metallicity gradients (as the gradient in the ISM is commonly believed to flatten out over time).

That is not to say that migration does not affect stars on hotter orbits, such as thick disk stars, as was shown in Solway et al. (2012), merely that stars on colder orbits may be preferred.

In Figure 9, we slice our data into four subsets on the basis of Z_max, find the corresponding metallicity gradients, and also plot their chemical and kinematic distributions, as in Figure 6. As we expected, the tighter Z_max constraints accentuate the trend for older stars to have flatter metallicity gradients.

Zoom In Zoom Out Reset image size

Two further trends become apparent here. First, we note that the chemical and orbital spaces are interlinked: the low-Z_max selected data in Figure 9 are also more α poor than the rest of the sample, while the high-Z_max data are more α rich. Second, we note that as we transition from the "thin disk" sample to the opposite "thick disk" sample, the gradient flattens out to become more similar at all ages and increases first toward zero, and then to positive values.

For intermediate-Z_max stars, we see a slight "uptick" with the youngest stars having a shallower [Fe/H] versus R gradient than intermediate-age stars. We believe this could be contamination from the thick disk (because the gradient is more positive than the most "thin disk" selection). While thick disk stars are not expected to be prominent in this age range, they could be "observationally young" yet "kinematically old" if they are indeed blue stragglers, as proposed by some authors and discussed above.

Note that a Z_max cut, particularly a strict one, could remove some of the strongest migrators in the data if they also flare as they migrate. This would mean that, for example, a high-metallicity star that migrates from the center to the outskirts, which should raise the [Fe/H] value in that outer bin, may be removed. This would work to reduce the metallicity in outer bins and increase the steepness of our gradients. Figures 9 and 10 can be investigated to consider the magnitude of this effect.

Zoom In Zoom Out Reset image size

Two populations may be seen in a bifurcation of the [α/Fe]–[Fe/H] space in the local neighborhood. In Figure 5 we showed the abundance plot for our data and noted a dense concentration of stars at low α abundances and a sparse cloud at high α abundances. We split these data into four sectional divisions. The strong concentration of stars with α abundances below 0.15 would be the conventional "chemical thin disk." Above this, at iron abundances below zero, would be the "chemical thick disk."

On the basis of the divisions shown in Figure 5, we again slice our data into four subsets and then fit the metallicity as a function of R_g in age bins in Figure 10.

We see that the α-poor stars inhabit the thin-disk-like areas in the orbital characteristic plots, with lower Z_max. They also have the strongest trend in [Fe/H] versus R_g as a function of age, similar to the orbitally selected thin disk. As we progress toward thick disk chemistries (α-rich stars), this gradient again moves toward positive values and flattens out.

For intermediate α-abundance stars, we again see an "uptick" where the youngest stars have a shallower gradient than intermediate-age stars.

An oddity in Figure 10 provides a nice opportunity to investigate how all of these figures fit together. One may note that, in the second column of Figure 10, the calculated gradients (panel 2) appear steeper for all ages of stars than the gradients shown in Figure 6 (panel 1), except the youngest bins, which are slightly shallower. Because the sample in this column is quite similar to the sample used in Figure 6, we wonder if there is a disagreement.

If we look to the first column in Figure 10 (panel 1), which is the second half of the chemical space covered by the sample in Figure 6, we can see that adding this sample to that of panel 2 will act to pull the gradient upwards at all points. If we consider the leftmost column (panel 1) of Figure 9, we can see that restricting our data to having a small Z_max (the other constraint in Figure 6) will act to pull the gradient for young stars downward.

In this way, Figures 6 and 7 can be reconstructed qualitatively from Figures 9 and 10.

As described, our determination of R_g is drawn directly from the rotation curve of our potential. As such, it is worth investigating the effects of different potentials on our results. We consider, in addition to the modified MWPotential2014 potential we have used for the main portion of this work, the potentials of Irrgang et al. (2013) and McMillan (2017) in Figure 11. The global trends in the fits of the metallicity gradient seem relatively unaffected by the choice of potential.

Zoom In Zoom Out Reset image size

While the choice of potential does not have much of an effect on the measured metallicity gradient, the choice of solar parameters—that is, the radius of the solar circle from the center of the Galaxy and local rotation speed—can make a relatively larger difference. In Figure 12, we show the effect of changing these parameters. The gradient will shift upward and downward in magnitude for various choices of the solar radius, R₀. The general trend and slope of the gradient in age in our thin disk sample remain relatively unchanged. The choice of rotation speed can change the trend more dramatically, altering the slope of the metallicity gradient with age. In all of these examples, however, the pattern of younger stars having steeper gradients than older stars remains.

Zoom In Zoom Out Reset image size

Age determinations are complicated. One particular weakness of the method that we use in this work, isochrone matching, is that stellar model isochrones tend to grow closer in the areas of the giant branch and the cooler portions of the main sequence. This causes isochrone age determinations to lose some of their discriminating power between ages in these areas, as the differences between models grow smaller compared to the typical stellar atmospheric parameter's errors.

Many works circumvent this issue by restricting their sample to the turnoff and subgiant branch areas of the color–magnitude diagram, because, in this area, isochrones differ more dramatically with age, making age determinations more reliable.

We therefore construct a smaller, high-quality sample from our data, consisting of 171,691 stars with 10% or better age errors as well as temperatures hotter than 5750 K, or log(g) < 4 and hotter than 5250 K. This should remove the populations most likely to have confused age estimates. The results of fitting the metallicity gradients of this smaller, more reliable sample are shown in Figure 13 and are similar to those found using our full sample.

Zoom In Zoom Out Reset image size

Star formation in the Milky Way is commonly thought of as growing in an inside-out (beginning in the central regions before expanding radially) and possibly upside-down (beginning with a vertically hot distribution that cools over time) fashion (for a more detailed explanation, see Bird et al. 2013). In this scenario, one would expect the oldest stars to have steep radial metallicity gradients (as the central regions are enriched and little to no star formation occurs in the periphery) and the youngest stars to have shallower gradients (as star formation and enrichment spreads outward through the extent of the modern-day disk).

Instead, we observe the opposite trend in our thin disk stars. This could be direct observational evidence of radial migration, the process by which stars' orbital radii may be shifted inward or outward in the galaxy and which results in a flattening of metallicity gradients over time.

Our analysis suggests that a single-parameter cut is not effective at selecting the population which gives the most constant flattening rate, because we see young-age flattening and "upticks," which lead to "U" shapes in the d[Fe/H]/dR_g versus age plots (the left and middle-left panels in Figures 9 and 10, as compared to Figure 6). This behavior could be indicative of contamination from other populations. Using kinematic and chemistry cuts in tandem seems to isolate the gradient-flattening trend more effectively than either individual cut.

If these upticks are caused by a contamination from a distinct population such as the thick disk (which seems to have a higher, near-constant gradient of 0.02 dex kpc⁻¹), it poses an interesting conundrum in that we see evidence of this contamination most in the youngest age bins. It should not be an effect of noisy observations "piling up" in the edge (youngest) age bins, because the trend extends past the edge of the isochrone grid to ages of 5 Gyr. That the chemical thick disk is not exclusively old is an idea that has drawn attention recently with growing numbers of observations of young, α-enhanced stars; however, many of these observations imply mass transfer origins, so their atmospheres could appear "young" even though they are members of a kinematically older population.

The selection that makes the "churning" signal the strongest—that is, the cut where the trend between d[Fe/H]/dR_g is most monotonic with age—is the cut based on orbital Z_max. This could be because of a "provenance bias," or preference for migratory forces tend to favor (although not exclusively) stars on vertically cooler and less eccentric orbits, as mentioned in Section 4.1.1.

The chemical and orbital thick disk selections yield positive metallicity gradients that are similar for all ages. That the gradients are similar, regardless of age, implies that the population is well mixed in age. It is surprising that it is not well mixed in radius, which would, presumably, result in a zero radial gradient. It is further surprising that the gradient is slightly positive. This might indicate that the chemical evolution of the thin and thick disks differs.

Positive gradients have been observed in external galaxies' thick disks (Cresci et al. 2010). In the Milky Way, this has been observed indirectly in the velocities of thick disk stars as a function of metallicity (as in Lee et al. 2011; Kordopatis et al. 2013, among others), where it is seen that more metal-rich thick disk stars actually rotate more quickly (implying that we are observing them near pericenter) than lower-metallicity thick disk stars.

Other studies, such as that by Nidever et al. (2014), find that the high-α stars in APOGEE have remarkably similar chemical distributions at all observational radii and heights from the plane, a result confirmed by Hayden et al. (2015). Loebman et al. (2016) posit that this effect could arise in a simulation as a result of these stars being formed in a narrow radial range in a limited time frame.

We suggest that these observations of no gradient in previous works studying high-α stars in R_present could be consistent with our observation of a slightly positive gradient in R_g for our "thick disk" selection if the data are blurred. It is worth mentioning that high-α stars are not necessarily synonymous with thick disk stars. We again note that selection functions make such comparisons more complex, especially when comparing the results of one spectroscopic survey to another.

It has been suggested that a positive d[Fe/H]/dR gradient may arise from the accretion of pristine gas into the center of the Galaxy at early epochs (Cresci et al. 2010; Curir et al. 2012) and that it may be a natural occurrence during inside-out formation (for example, Schönrich & McMillan 2017; Kawata et al. 2018).

Our findings therefore are in agreement with other observations with respect to the positive gradients (for example Jia et al. 2018, see also Table 1 of Önal Taş et al. (2016)). That this gradient is similar at all ages is a result that should help constrain the formation history of the Milky Way's thick disk.