THE SPATIAL STRUCTURE OF MONO-ABUNDANCE SUB-POPULATIONS OF THE MILKY WAY DISK

Fitting the spatial-density profiles of various G-dwarf sub-samples must account for the fact that the observed star counts do not reflect the underlying stellar distribution, but are strongly shaped by (1) the strongly position-dependent selection fraction of stars with spectra (see Figure 11), (2) the need to use photometric distances that in turn depend on the color and metallicity distribution of the sample (as the magnitude-limited SEGUE sample corresponds to a color- and metallicity-dependent distance-limited sample), and (3) the pencil-beam nature of the SEGUE survey. To properly take all of these effects into account, we need to use forward modeling: in what follows we fit stellar-density models to the data by generating the expected observed distribution of stars in the spectroscopic sample, based on our model for the SEGUE selection function and the photometric-distance relation; this predicted distribution is then compared to the observed star counts. We show below how this can be expressed as a maximum likelihood problem. This general density-fitting methodology applies to any spectroscopic survey, with minor modifications, and needs to be applied to obtain selection-corrected distributions from spectroscopically selected stellar samples. In particular, this methodology needs to be applied to constrain the structural parameters of abundance-selected samples in the Milky Way.

As the photometric distance estimates depend on the g − r color, metallicity [Fe/H], and apparent r-band magnitude, and because the selection function is a function of position, r, and g − r, we need to model the observed density of stars in color–magnitude–metallicity–position space, λ(l, b, d, r, g − r, [Fe/H]). This density of stars can be written as

$$\begin{eqnarray} &&\lambda (l,b,d,r,g-r,[\mathrm{Fe/H}])= \rho (r,g-r,[\mathrm{Fe/H}]|R,Z,\phi)\nonumber\\ &&\qquad\times \nu _*(R,Z,\phi) \times |J(R,Z,\phi;l,b,d)|\nonumber\\ &&\qquad \times S(\mathrm{plate},r,g-r).\quad \end{eqnarray} \tag{ 4 }$$

Here, (R, Z, ϕ) are Galactocentric cylindrical coordinates corresponding to rectangular coordinates (X, Y, Z), which can be calculated from (l, b, d). The factor ρ(r, g − r, [Fe/H]|R, Z, ϕ) is the number density in magnitude–color–metallicity space as a function of position (see further discussion in Appendix B). The |J(R, z; l, b, d)| is a Jacobian term because of the (X, Y, Z) → (l, b, d) coordinate transformation; the crucial factor S(plate, r, g − r) is the selection function as given in Equation (A2). Finally, ν_*(R, Z, ϕ) is the underlying spatial number density of the sample; we stress that this is a density as a function of rectangular coordinates (X, Y, Z) that we evaluate through (R, Z, ϕ), i.e., its dimension is 1/(spatial unit)³. In what follows we will assume that our models for this density (e.g., exponentials in the vertical and radial directions) are characterized by a set of parameters denoted as θ and that the density is axisymmetric, such that ν_* ≡ ν_*(R, Z|θ).

The likelihood of a given model for the density ν_*(R, z|θ) is given by that of a Poisson process with rate parameter λ,

$$\begin{eqnarray} \ln \mathcal {L} &=& \sum _i [ \ln \lambda (\lbrace l,b,d,r,g-r,[\mathrm{Fe/H}]\rbrace _i|\theta) ]\nonumber\\ && -\int {d}l \,{d}b\, {d}d\, {d}r\, {d}(g-r) \,{d}[\mathrm{Fe/H}]\nonumber\\ &&\times \lambda (l,b,d,r,g-r,[\mathrm{Fe/H}]|\theta), \end{eqnarray} \tag{ 5 }$$

where the integral is over the domain surveyed and i indexes the observed objects. Because the Jacobian, the selection function, and the density in magnitude–color–metallicity space only enter λ multiplicatively (Equation (4)), their contribution to the first term (ln λ) in Equation (5) is a constant that does not depend on the density parameters. Thus, up to a term that does not depend on θ, the log likelihood is equivalent to

$$\begin{eqnarray} \ln \mathcal {L} &=& \sum _i [ \ln \nu _*(R,z|\lbrace l,b,d\rbrace _i,\theta) ]-\int {d}l \,{d}b\, {d}d\, {d}r\, {d}(g-r) \,\nonumber\\ &&\times {d}[\mathrm{Fe/H}]\,\lambda (l,b,d,r,g-r,[\mathrm{Fe/H}]|\theta) \,. \end{eqnarray} \tag{ 6 }$$

Note that the Jacobian, the density in the magnitude–color–metallicity space, and the selection function only enter through the second term, and do not need to be evaluated on a star-by-star basis. The second term in Equation (6)—the normalization integral—can be written as (assuming that the density does not depend on (l, b) over the area of a plate, although this can easily be relaxed)

$$\begin{eqnarray} &&\int {d}l\, {d}b\,{d}d\, {d}r\, {d}(g-r) \,{d}[\mathrm{Fe/H}]\,\lambda (l,b,d,r,g-r,[\mathrm{Fe/H}]|\theta) \nonumber\\ &&= A_p\,\sum _{\mathrm{plates}\ p} \int {d}(g-r)\,{d}[\mathrm{Fe/H}]\,{d}r \, S(p,r,g-r) \nonumber\\ && \times\int dd\, \rho (r,g-r,[\mathrm{Fe/H}]|R,Z,\phi)\,d^2\,\nu _*(R,z|l,b,d,\theta)\,, \end{eqnarray} \tag{ 7 }$$

where A_p is the area of a SEGUE plate (approximately 7 deg²).

In the following, we analytically marginalize over the amplitude of the rate λ with a logarithmically flat prior. In that case, the log likelihood becomes

$$\begin{eqnarray} \ln \mathcal {L} &=& \sum _i [ \ln \nu _*(R,z|\lbrace l,b,d\rbrace _i,\theta) -\ln \int {d}l \,{d}b\, {d}d\, {d}r\, \nonumber\\ &&\times{d} (g-r) \,{d}[\mathrm{Fe/H}]\,\lambda (l,b,d,r,g-r,[\mathrm{Fe/H}]|\theta)]. \quad\ \ \ \end{eqnarray} \tag{ 8 }$$

Note that the normalization integral is now moved inside of the logarithm.

In Appendix B, we discuss how we include the magnitude–color–metallicity factor ρ(r, g − r, [Fe/H]|R, Z, ϕ) in the likelihood.

We fit number-density models for the various abundance sub-populations, consisting of a disk with an exponential profile in both the vertical and radial directions, plus a constant density

$$\begin{eqnarray} \nu _*(R,Z) &=& N(R_0)\bigg[\frac{1}{2\,h_z}\, \exp \left(-\frac{R-R_0}{h_R}\right)\,\nonumber\\ &&\times \exp \left(-\frac{|Z|}{h_z}\right) + \frac{\beta _c}{24}\bigg]\,, \end{eqnarray} \tag{ 9 }$$

where N(R₀) is the vertically integrated number density at R₀. We refer to this model below as a single-exponential disk fit, as in all cases the data imply β_c ≪ 1. We also fit combinations of exponential disks as

$$\begin{eqnarray} \nu _* (R,Z) &=& N(R_0)\,\left[\frac{1-\beta _2}{2\,h_z}\,\exp \left(-\frac{R-R_0}{h_R}\right)\,\exp \left(-\frac{|Z|}{h_z}\right)\right. \nonumber\\ &&\left.+\, \frac{\beta _2}{2\,h_{z,2}}\,\exp \left(-\frac{R-R_0}{h_{R,2}}\right)\,\exp \left(-\frac{|Z|}{h_{z,2}}\right)\right]\,. \end{eqnarray} \tag{ 10 }$$

In particular, in the Z-direction, this is analogous to traditional density fits based on photometric data, which require (at least) two exponential components. We do not fit for the overall normalization, N(R₀), as we are interested primarily in the shape of the stellar-density profile.

To determine the best-fit parameters and their uncertainties we use Powell's method for minimization (Press et al. 2007), and then MCMC-sample the posterior distribution function, obtained by multiplying the likelihood in Equation (8) with flat logarithmic priors for the scale parameters (h_z, h_R, h_{z, 2}, h_{R, 2}) and flat priors on the contamination-fraction parameters (β_c, β₂), using an ensemble MCMC sampler (Goodman & Weare 2010; Foreman-Mackey et al. 2012).

In Appendix D, we discuss tests of the fitting methodology on mock data sets made up of single-exponential disk components observed using the SEGUE sampling. These tests show that we can recover the input density structure to within the MCMC-determined uncertainties over the range of inferred scale heights, scale lengths, and sample sizes found below.

For the α-old sample, the fit results for single-exponential profiles in R and Z, and for a combination of two exponential profiles for both R and Z, are given in Table 1. The model with two exponentials in both R and Z is preferred, but the parameters of the dominant double-exponential disk are similar for both fits. That is, even when we give the model the additional freedom of two vertical scale heights, the data lead us to employ only a single-exponential scale height. There is no evidence for a thinner component in the α-old abundance range. We see that the α-old sample is dominated by a population of stars with a scale height of 686 ± 11 pc, and a short scale length of 2.01 ± 0.05 kpc (consistent with the rough estimate of 2 kpc based on a handful of stars by Bensby et al. 2011 and the indirect dynamical estimate of 2.2 ± 0.35 kpc of Carollo et al. 2010).

We have split the α-old sample into more metal-poor and more metal-rich sub-samples by cutting the sample at [Fe/H] = −0.7. This is close to the median [Fe/H] of the α-old sample. The metal-poor sub-sample may be identified with the metal-weak thick-disk (MWTD) population discussed by Carollo et al. (2010), which they argue covers the metallicity range −1.8 ⩽ [Fe/H] ⩽ −0.8. The resulting fits for the spatial structure of these sub-samples are given in Table 1. The inferred scale lengths for these sub-samples are equal to within the uncertainties. However, the scale height of the more metal-poor sample is 856 ± 20 pc while that of the more metal-rich sample is 583 ± 16 pc. The radial scale length of the MWTD determined from the indirect dynamical analysis of Carollo et al. (2010) is roughly 2 kpc, while the scale height is 1.36 ± 0.13 kpc.

We have also split the α-old sample into two bins in [α/Fe] by splitting the sample at [α/Fe] = 0.35. The best-fit density profiles, given at the bottom of Table 1, again have similar scale lengths, around 2 kpc, and different scale heights. The stars that are most enhanced in α-elements have the largest scale height (765 ± 15 pc) and the shortest scale length (1.89 ± 0.04 kpc), while the less α-enhanced stars have a smaller scale height (627 ± 18 pc) and longer scale length (2.23 ± 0.1 kpc). As the latter dominate the full α-old sample, their scale height is very similar to that inferred for the full sample. We explore the dependence of the disk parameters on [Fe/H] and [α/Fe] in more detail in Section 4.3 below.

The results for single-exponential disk fits and double-exponential disk fits for the α-young sample are given in Table 2. The double-exponential disk fit model is formally preferred, but the parameters of the dominant double-exponential disk are again similar for both fits. We see that the α-young sample is dominated by a population of stars with a low scale height of 256 ± 4 pc and a long scale length of 3.6 ± 0.2 kpc.

The second double-exponential disk in the best-fit model for the α-young sample has a scale height of 664 ± 132 pc, which is consistent with the scale-height measurement of the α-old sample above. However, the fraction of stars in this secondary component is too small to constrain its scale length, and is conceivably simply a result of "abundance contamination" of the sample.

Density fits for α-young samples with the same [α/Fe] limits as the nominal α-young sample shown in the top panel, but that are more metal-poor, are also given in Table 2. We do not measure any radial density decline for these more metal-poor α-young samples, and short scale lengths for these samples are ruled out by the data. We consider this further in Section 4.3 and in the discussion section below.

When we split the α-young sample into two pieces, by cutting at [α/Fe] = 0.15, we find that the more α-enhanced sample has the shortest scale length (2.3 ± 0.2 kpc) and the largest scale height (348 ± 13 pc). The sample with [α/Fe] closer to solar has a longer scale length of 4.3 ± 0.2 kpc and a smaller scale height of 239 ± 4 pc.

In the previous two sections, we found that sub-samples of stars defined by their element abundances appear to have a simple spatial structure, approximated by a single exponential in the radial and vertical directions. The scale lengths and heights of these sub-sets seem to vary systematically with the abundances: the α-old sample has a shorter scale length than the α-young sample, and if we split those two samples further in [α/Fe], the part of the α-young sample that has the closest to the solar [α/Fe] ratio has the longest scale length and the smallest scale height. We also noticed that populations with [α/Fe] <0.25 have longer scale lengths and scale heights with decreasing [Fe/H].

To further investigate these trends, we have fitted disk models with single-exponential profiles in R and Z to sub-populations of stars with narrow bins in [Fe/H] and [α/Fe]. We divide stars into bins of width 0.1 dex in [Fe/H] and 0.05 dex in [α/Fe], and only fit those bins with more than 100 stars. The results from these fits are shown in Figure 4. The populations in the lower left part of the [α/Fe]–[Fe/H] diagram all have best-fit scale lengths in excess of 4.5 kpc.

Zoom In Zoom Out Reset image size

We also fitted two-component models, i.e., two exponential disks, to each of the bins, but found that these led to overfitting, and only marginal improvements in the likelihood for the best fit. Thus, for narrow bins in elemental-abundance space, the sub-populations are very well described by single-exponential profiles in the R- and Z-directions.

A different view of the results in Figure 4 is given in Figure 5. The results in the different [Fe/H]–[α/Fe] bins are shown as a function of scale length and scale height; the points are color-coded by their [α/Fe] or [Fe/H] dependence, and the size of the points corresponds to the total stellar surface-mass density—corrected for mass and sample selection effects—in each population (calculated in Bovy et al. 2012a). Figure 5 also shows the uncertainties in the inferred parameters; the formal uncertainty in the scale height for some points is so small that it cannot be seen. The bins with dashed error bars lie in a part of the abundance plane where abundance contamination is likely to be the most severe, where the [α/Fe]-based age ranking is least reliable, and where the spatial properties change most rapidly. They contain <5% of the disk surface mass.

Zoom In Zoom Out Reset image size

We see that these fits for mono-abundance sub-components flesh out the main trends we noted in the broader [Fe/H] and [α/Fe] ranges above. At any given metallicity [Fe/H], the scale length increases and the scale height decreases when moving from α-old to α-young populations. At any given α-age, the scale length and the scale height increase for the more metal-poor components, implying an outward metallicity gradient. And, as Figure 5 shows most clearly, increasing scale lengths are correlated with decreasing scale heights (except for a few bins on the boundary between the very long scale lengths at low metallicity and solar α-enhancement and the shorter scale lengths of the α-old populations; see further discussion in Section 5.5). From Figure 4 it is clear that neither [α/Fe] nor [Fe/H], on its own, accounts for the trends in scale height and scale length. We discuss what this implies for disk formation and evolution in Sections 5.4 and 5.5, respectively.

Figure 6 shows the results of fitting two components with exponential profiles in both R and Z to each abundance bin. The scale height of the dominant component is shown against the best-fit scale height, when fitting a single-exponential profile in R and Z. We see that these scale heights are strongly clustered around the one-to-one correspondence line. Thus, for each bin, a single vertical exponential suffices to explain the observed number counts. The fact that the two measurements agree better than would be expected, given the uncertainties shown, is due to the fact that the scale heights for each bin are strongly correlated when fitting a single- or a double-exponential profile in R and Z. Overall, Figure 6 confirms that a single-exponential model in Z and R is a good model for the spatial structure of mono-abundance sub-populations.

Zoom In Zoom Out Reset image size

In Appendix D, we perform a test to determine whether abundance uncertainties can plausibly lead us to find spurious disk components between a "thin" and a "thick" component. That is, we ask whether it is plausible that an underlying density dominated by distinct thin- and thick-disk components can be smoothed by abundance errors into the density structure we inferred in Figures 4–6. This test shows that if this were the case, every bin is preferentially fitted with two components, corresponding to the input thin and thick components. The equivalent of Figure 6, shown in the bottom right panel of Figure 19, is qualitatively different, with a distinct difference between the single-component scale height and the scale height of the dominant component in the two-component fit.

To test whether the analysis in this section is influenced by our signal-to-noise ratio cut of S/N >15, we have repeated the analysis with a cut of S/N >30, as also used by Lee et al. (2011b). The equivalents of Figures 4–6 look qualitatively the same, albeit with larger uncertainties for each bin, and the dependence of h_z and h_R on elemental abundance is the same as that inferred from the sample with the S/N >15 cut. The number of ([Fe/H], [α/Fe]) bins with more than 100 stars is smaller, but the inferred (h_z, h_R) for those bins with more than 100 stars when using the S/N >30 cut are consistent within the uncertainties with those found with the less restrictive S/N cut. We stress that even when selecting stars with S/N >30, the equivalent of Figure 6 does not show any sign of a second component in the mono-abundance bins.

To perform the binning in this section, we used narrow bins of 0.1 dex in [Fe/H] and 0.05 dex in [α/Fe]. These bins are somewhat narrower than the total typical uncertainty (≈0.15 dex in [Fe/H], ≈0.07 dex in [α/Fe]; Bovy et al. 2012b), but we prefer to oversample, rather than undersample, to avoid smoothing out underlying structure. The analysis in each bin holds irrespective of the bin size. What matters for the analysis is that the data in each bin are disjoint, such that the bins are statistically independent.

The absolute values of the distance scales measured in this paper are subject to distance systematics, which we discuss in this subsection. We have used the data-driven photometric-distance relation from Ivezić et al. (2008) to infer the spatial structure of the various samples of stars, but an alternative photometric-distance relation can be obtained by using the An et al. (2009) stellar isochrones in the SDSS passbands. These isochrones depend on [Fe/H] as well as on [α/Fe], although in practice a linear relation between [α/Fe] and [Fe/H] is assumed, and the spectroscopically measured [α/Fe] is not used directly to estimate the photometric distance. In the top panel of Figure 1, we compared the distance moduli derived using the An et al. (2009) stellar isochrones with those obtained using the Ivezić et al. (2008) relation for a few values of [Fe/H]. We see that, for the values of [Fe/H] that span most of our sample, the distance modulus difference is −0.2 mag, corresponding to a systematic difference in the inferred distances of about 9%, nearly independent of color. Thus, if we had used the An et al. (2009) photometric distances, we would have obtained scale lengths and scale heights that were 9% shorter.

A second distance systematic that could influence our results is the Malmquist bias (Malmquist 1920, 1922)—the fact that brighter stars are overrepresented in a magnitude-limited survey. For our relatively bright sample, this is dominated by the finite width of the photometric-distance relation. The Malmquist bias in absolute magnitude is apparent-magnitude dependent and approximately equal to −σ²dln A(r)/dr, where A(r) is the differential number count as a function of apparent magnitude and σ is the dispersion in the absolute magnitudes (either due to photometric uncertainties or due to intrinsic scatter in the photometric-distance relation). Conservatively assuming that the combination of the finite width of the photometric-distance relation and the photometric uncertainties is 0.2 mag, and that the underlying density is constant, the Malmquist bias would be of order 2.5%. However, due to the exponential falloff of the density in both the R- and Z-directions, the differential number counts are (1) flat near the peak induced by the vertical exponential and (2) for most apparent magnitudes |dln A(r)/dr| is less than 1. Therefore, the Malmquist bias is at most about 2%, and will not strongly affect the measurement of the vertical scale height in particular.

We have assumed throughout our analysis that all of the stars in our sample are single. The presence of unresolved binaries will lead us to underestimate scales, as these binaries will appear to us as brighter, and thus closer, single stars. The binary fraction and companion-mass distribution for G-type dwarfs remain controversial, but it appears that the overall binary fraction for G dwarfs is approximately 40% (Abt & Levy 1976; Duquennoy & Mayor 1991; Raghavan et al. 2010), similar to but slightly larger than that of M dwarfs (Fischer & Marcy 1992; Raghavan et al. 2010). The distribution of companion masses is poorly known, and could range from being peaked around 20% of the primary's mass (Duquennoy & Mayor 1991), to being relatively flat between 20% and 100% of the primary's mass (Raghavan et al. 2010), with numerical simulations indicating that multiple-star systems form preferentially with approximately equal-mass members (Bate 2005), and an overall multiplicity fraction of around 40% (Bate et al. 2003). Lower-metallicity stars most likely have a higher binary fraction (Machida et al. 2009), and could reach 100% for [Fe/H] <−0.8 (Raghavan et al. 2010).

For a likely scenario where 40% of our α-young sample is made up of binary stars (ignoring higher-order multiplicities) with a flat distribution of companion masses between 20%  and 100% of the primary's mass, the magnitude would be overestimated on average by 0.12 mag, such that the scale height and scale length would be underestimated by about 6%. If 70% of the α-old sample were to consist of binary systems (taking into account the rising binary fraction with decreasing metallicity), the magnitudes would be overestimated by approximately 0.21 mag, and the α-old scale heights and scale lengths would be underestimated by 10%. These biases are somewhat larger than the statistical uncertainties on our results, but they are similar to the overall distance-scale uncertainty (see above), and they do not change the conclusion that the α-old scale length is much shorter than that of the α-young sample. Even in a worst-case scenario, where all binary systems have equal-mass companions and where 100% of the α-old stars are in binaries, the α-old scale length would still be ≲ 2.8 kpc (40% up from 2 kpc), which is shorter than the scale length measured for the α-young sample in Table 2 and Figure 5, and the α-young scale lengths themselves would also increase by about 15% in this scenario. In principle, a careful spectral analysis of the SEGUE spectra themselves could provide direct constraints on the (unresolved) binary contamination in this sample.

In our density fits we have mostly fitted disk components to the data, except for the single-exponential disk model where we added a uniform density (Equation (9)). We thus assumed that the stellar halo does not influence our disk fits, beyond what can be described by a uniform density across our survey volume. We can estimate the expected number of halo stars in our sample using the Bell et al. (2008) density fits to the smooth stellar halo. We run the Bell et al. (2008) stellar-halo density through the G-star SEGUE selection function, and marginalize over g − r color using a flat distribution over 0.48 ⩽ g − r ⩽ 0.55, and over [Fe/H] using the Ivezić et al. (2008) halo metallicity distribution (mean [Fe/H] = −1.52, width = 0.32). We then find that for ≈10⁸ G-type stars between 1 and 40 Galactocentric kpc in the stellar halo, there should be about 100 halo stars in our sample, compared to the total sample size of 30,353 G-type dwarfs. Hence, the halo contamination is very small and does not influence the fits. Additionally, halo contamination will be most severe for the α-old sub-populations, and this contamination should work to increase the inferred scales (length and height). Therefore, the result that the radial scale length of α-old sub-populations is shorter than that of α-young sub-populations is robust against any halo contamination.

The density fits in this paper are the first to constrain the vertical scale height and radial scale length of numerous disk sub-components, defined using elemental abundances alone, from a large sample of stars. Our results show that the vertically thicker-disk sub-components—when chemically defined—have a much shorter scale length than the thinner-disk sub-components, which is opposite to traditional purely geometric disk decompositions (e.g., Robin et al. 1996; Ojha 2001; Larsen & Humphreys 2003), which typically find that the thick-disk component has a longer scale length than the thin disk, and that this scale length is ≳ 3.5 kpc (e.g., Jurić et al. 2008).

When we fit the spatial structure in our approach, taking stars of all metallicities (specifically, the combination of our α-old ("thick") and α-young ("thin") samples), we can recover the result of purely geometric decompositions: the thin-disk component—i.e., the component with the lowest scale height, ≈300 pc—gets paired with the shortest scale length (≈2 kpc), while the thicker-disk component gets assigned both the largest scale height and scale length (for our particular sample fit with a combination of three double-exponential disks these are ≈600 pc and ≈2.4 kpc, with a small component with an even larger scale height and scale length). Thus, it seems that purely geometric decompositions naturally associate the longest scale length with the largest scale height. That both geometrically determined scale lengths are shorter than the scale length of the α-young sample is due to the fact that the metallicity distribution for the entire sample extends down to [Fe/H] = −1.5, such that the model "expects" many low-metallicity stars in the "thin" component at large distances (as the model does not contain the information that the thin component has higher metallicities), which are not observed. Therefore, metallicity and α-element abundances, which are manifestly quantities that can identify sub-samples of stars independent of their spatial structure and kinematics, lead to a qualitatively different decomposition into two (or more) sub-components than the purely geometrical approach, with its inherent risk of circular reasoning.

The distinctive changes of the global disk structure with abundance, especially with the age proxy [α/Fe], should provide valuable clues to the formation of the Milky Way's disk. While a concrete and quantitative model comparison is beyond the scope of this paper, we discuss some of the qualitative implications here. As mentioned in Section 1, the overall radial-density profile of the stellar disk is expected to be conserved even in the face of large-scale radial migration, but the radial profile of sub-components will tend to relax to the mass-weighted mean radial profile. Thus, a difference in the radial distribution of various populations of stars today is a less-pronounced version of more different initial radial distributions (at formation). Assuming that the [α/Fe] ratio is an adequate proxy for age (e.g., Schönrich & Binney 2009b), our results then imply that the α-enhanced, hence oldest, populations are more centrally concentrated—have a shorter scale length—than populations with α-abundances that are closer to solar, and therefore younger. This is direct observational evidence for inside-out formation of galactic disks across the presumed age range of our sample, 1–10 Gyr, where the inner parts of the disk form before the outer part of the disk. A similar age dependence of the exponential scale length has been found in several external galaxies (de Jong et al. 2007; Radburn-Smith et al. 2012).

Second, our analysis shows that our Milky Way has not only a metallicity gradient among its youngest stars, but that it has always had one (Cheng et al. 2012): at a given [α/Fe], standing in for age, sub-populations with lower [Fe/H] have a longer scale length than more metal-rich stars. This picture is confirmed by looking at the orbital properties of the stars when integrating our sample of G-type dwarfs in a simple model for the Milky Way's potential, made up of a Miyamoto–Nagai disk with a radial scale of 4 kpc and vertical scale of 300 pc contributing 60% of the radial force at the solar radius, a Hernquist bulge with a scale radius of 600 pc contributing 5% of the radial force, and a Navarro–Frenk–White halo with a scale radius of 36 kpc that contributes 35% of the rotational support at the Sun's position. The median of the mean orbital radii as a function of elemental abundance is shown in Figure 7 (see also Lee et al. 2011b; Liu & van de Ven 2012 for similar figures of the eccentricity and rotational velocity). We see that stars with [α/Fe] <0.25 and lower [Fe/H] are thin-disk stars that live, on average, farther out than more metal-rich stars. Thus, the longer scale length for outer-disk stars, combined with the fact that, for solar [α/Fe], decreasing [Fe/H] is correlated with decreasing age (e.g., Schönrich & Binney 2009b), implies that the outer part of the disk formed later than the inner part.

Zoom In Zoom Out Reset image size

We have assumed that [α/Fe] is an adequate proxy for age, such that the mono-abundance populations that are more [α/Fe]-enhanced are older than the populations with solar [α/Fe]. This is typically the case in standard scenarios for the star formation history of the Milky Way disk, in which [α/Fe] steeply drops around 2–3 Gyr due to the onset of Type Ia supernovae (Dahlen et al. 2008; Maoz et al. 2011), and then stays roughly constant, although the value of [Fe/H] at which the [α/Fe] downturn happens depends on the star formation history (Matteucci & Recchi 2001). Only if the local star formation was characterized by bursts of star formation can younger populations of stars have similar levels of [α/Fe] as older stars (Gilmore & Wyse 1991). Most current fits of the local star formation history prefer a smooth history (e.g., Aumer & Binney 2009), although it is difficult to rule out epochs of enhanced star formation (e.g., Rocha-Pinto et al. 2000).

The spatial structure inferred for mono-abundance sub-populations (Figures 4 and 5) shows two important results: first, there is a tight anti-correlation between the scale heights and scale lengths of the sub-components. Secondly, there is a continuous distribution in scale height when moving from α-enhanced, metal-poor populations to stars with solar [α/Fe] and [Fe/H]. This suggests that the α-old, "thick,"and the α-younger, "thin," regime of the stellar disk are not two separate entities, but merely opposite ends of the disk evolution spectrum (suggested before in Norris 1987, but never directly measured as we do here). This issue, which requires proper stellar-mass weighting of the sub-components, is worked out in a separate paper (Bovy et al. 2012a). Taken together, these findings suggest a continuous evolutionary mechanism created the observed scale-height distribution, rather than a discrete external heating or accretion event. Radial migration is an obvious candidate for this internal evolution mechanism. That the most centrally concentrated component of the disk is not only the (α-)oldest part, but also has the largest scale height, is a nearly inevitable condition, and hence a natural prediction, of any scenario where much of the disk scale-height distribution is created through radial migration. The α-old sub-population not only had the most time to evolve, but its centrally concentrated parent population implies that stars at 6 kpc < R < 12 kpc have migrated out by the largest factor.

A different internal explanation for the thicker-disk components in the Milky Way is that, rather than being thickened over the history of the Galactic disk, thick components were created thick during an early, turbulent phase in the formation of the disk (e.g., Bournaud et al. 2009; Förster Schreiber et al. 2009). If such a scenario is combined with a inside-out growth of the disk, and the disk remains turbulent over a significant fraction of its history, this formation scenario could plausibly explain the continuous dependence of disk structure on elemental abundance found in this paper.

Our result that the transition between the α-young, "thin," components and the α-old, "thick," components is smooth, rather than showing a clear separation between thin and thick components, may appear to be in conflict with local, high-resolution spectroscopic samples of stars (e.g., Reddy et al. 2006; Fuhrmann 2011; Navarro et al. 2011) or other analyses of the SEGUE data (e.g., Lee et al. 2011b). A detailed comparison between these and our results requires careful accounting for the spectroscopic volume sampling, which has not been done in the Lee et al. (2011b) analysis or for the high-resolution samples, except for the sample of Fuhrmann (2011), which is volume complete out to 25 pc. Without taking the volume selection into account, the sample used here also displays a bi-modality in the [Fe/H]–[α/Fe] plane (see Figure 2). We discuss this issue in more detail in Bovy et al. (2012a), but we note here that the apparent bi-modality in the observed number density of stars disappears when properly correcting for the spectroscopic sampling. Furthermore, the local, high-resolution analyses cannot directly measure the spatial distribution of stars of different elemental abundances (e.g., Fuhrmann 2011, which only has 15 high-[α/Fe] stars out to 25 pc; Reddy et al. 2006) and therefore rely on kinematics to argue that the vertical distribution of stars in the solar neighborhood is characterized by a bi-modal "thin"–"thick" disk dichotomy. This interpretation is driven by the selection of stars that are disjoint in [α/Fe], which leads to disjoint kinematics because the kinematics is a strong—and smooth—function of abundance as well (Bovy et al. 2012b). While the stellar content of different survey volumes can (and should) be connected by dynamics, we note that the effective volumes sampled by, e.g., the Fuhrmann (2011) survey and by our analysis differ by a factor of about 10⁵; hence, the extrapolation from one to the other is enormous. The analysis of the vertical kinematics of stars in our sample confirms the existence of the intermediate populations with scale heights between 400 and 600 pc and vertical-velocity dispersions of 30–35 km s⁻¹ (Bovy et al. 2012b).

Our finding that the scale length does not behave as smoothly as the scale height, as a function of [α/Fe], is presumably a consequence of the disk's formation history: here the increasing metallicity as a function of time (i.e., youth) and the radial metallicity gradient compete. As the mapping between [α/Fe] and age is not linear, but rather [α/Fe] steeply drops around 2–3 Gyr due to the onset of Type Ia supernovae (Matteucci & Recchi 2001; Dahlen et al. 2008; Maoz et al. 2011) and then stays roughly constant, the scale length should change similarly rapidly with [α/Fe]. The scale height, however, is determined by subsequent evolution, where radial migration transports stars to larger Galactocentric radii, where the lower disk density allows them to travel farther from the plane. Since this evolution is continuous, rather than sudden, and includes additional contributions from heating, trends in scale height versus elemental abundance should be expected to be smoother, even if radial migration is not the disk's dominant evolutionary mode. Our results are therefore consistent with a scenario where the thick-disk component is the inner part of the disk that formed at the earliest time, and either by having formed thick or through the effect of radial migration, has a large scale height at the present time.

A gas-rich merger, followed by intense star formation at an early time, could have affected the formation of the early disk (Brook et al. 2004), as seems consistent with the observed distribution of eccentricities of the thick-disk component (Sales et al. 2009; Dierickx et al. 2010; Wilson et al. 2011). However, it would lead to a scale length for the thicker component that is larger than that of the thinner component (Qu et al. 2011). It is clear that internal mechanisms must have played an important role during the evolution of the disk. However, we caution that the radial and vertical consequences of neither radial migration nor turbulent disk evolution, nor of satellite thickening, have been worked out in quantitative detail, and, in particular, resonant coupling between satellites and the disk might induce some similar observational signatures to radial migration.

The rapid change in the mean stellar population in an [α/Fe]–[Fe/H] abundance bin at the onset of Type Ia supernovae is also likely the explanation for the presence of the few points of intermediate [α/Fe] and [Fe/H] in Figure 5 that do not follow the anti-correlation between scale height and scale length; these bins, which do not contribute significantly to the total stellar mass (indicated by the size of the symbols in Figure 5), are also the bins that fall short of the one-to-one correlation between single- and two-disk fits in Figure 6. This provides further evidence of the fact that at the rapid [α/Fe] (age) transition our bins do not adequately resolve single components.

To determine the spatial distribution of the G dwarfs, we require a good understanding of the SEGUE G-star selection function, i.e., the fraction of stars that has been targeted by SEGUE and produced good enough spectra to derive the parameters needed in the present (or any other) analysis (e.g., S/N >15), and we need this selection fraction as a function of position, color, and apparent magnitude. The observed density of G-type stars is simply the product of the underlying density with the sampling selection function, suggesting that one constrains this underlying density by forward modeling of the observations.

The spectroscopic G-star target type was selected uniformly from the set of objects in the G-star color–magnitude box in the area and apparent magnitude range of the spectroscopic plug-plates (simply "plates" hereafter); thus the selection function can be reconstructed. The SEGUE survey implementation distinguishes between "bright" and "faint" plates, with bright plates containing stars with r ⩽ 17.8 mag and faint plates containing stars with r > 17.8 mag. For the purposes of the selection function, we assume that this separation at 17.8 mag is a hard cut, even though in reality some stars were observed on both bright and faint plates for calibration purposes, and some "bright" stars are part of faint plates, and vice versa, because of changes between the photometry used for target selection and that released as part of the SDSS DR7, which we employ here. Duplicates are resolved in favor of the higher S/N observation (typically on the faint plate as this has a longer integration time). We retain stars with r ⩾ 17.8 mag when they were observed on a bright plate, and we keep objects with r < 17.8 mag when they were observed on a faint plate, even though this should not happen in our model for the SEGUE selection function below. A total of 586 stars in the α-old sample and 47 stars in the α-young sample fall into this category; they do not influence any of the fits or conclusions in this paper.

We select the superset of targets by querying the SDSS DR7 imaging CAS⁸ for all potential targets in the color–magnitude box of the G-star target type in the area of a SEGUE plate (Yanny et al. 2009). These objects are primary⁹ detections (removing duplicates and objects from overlapping imaging scans) with stellar point-spread functions (PSFs) (type equal to 6). Objects must not be saturated, nor be close to the edge, nor have an interpolated PSF (interp_psf), and must not have an inconsistent flux count (badcounts). Furthermore, if the center is interpolated (interp_center), there should not be a cosmic ray indicated (cr). See Stoughton et al. (2002) for a description of the SDSS photometric flags. Using the superset of targets we determine for each plate the fraction of stars that were observed spectroscopically of all available targets.

To infer the dependence on color and apparent magnitude of the selection function, we look at the distribution of the potential G-star targets in color–magnitude space. This is shown in Figure 8. The distribution of the spectroscopic sample is overlaid. This shows that the spectroscopic sampling is relatively fair in g − r color, with some frayed edges because of changes between target and current photometry, and that the selection as a function of r-band magnitude tapers at the faint end, as should be expected when using an S/N cut. If all SEGUE plates were integrated to the same depth, the S/N cut should be a clean cut in r, but it is clear from Figure 8 that this is not the case. To distinguish between relatively shallow and relatively deep plates, we introduce the overall plate signal-to-noise ratio plateSN_r,

$$\begin{equation} \tt{plateSN\_r} = (\mathrm{sn1\_1}+\mathrm{sn2\_1})/2\,, \end{equation} \tag{ A1 }$$

where sn1_1 and sn2_1 are the r-band plate S/N for the two SDSS spectrographs (see Table 17 in Stoughton et al. 2002). The faintest spectroscopic G-type star per plate as a function of plateSN_r for the faint plates is shown in Figure 9 for the faint plates. This figure shows that there is a clear difference in the faintest object that could have been successfully observed at S/N > 15 between relatively shallow and relatively deep plates. The bottom panel of Figure 9 shows the S/N of stars on four plates chosen to cover a range in the overall plate S/N. This shows that the S/N > 15 cut for the entire sample translates into a fairly sharp r-band cut for each individual plate.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Our model for the SEGUE G-star selection function is then the following: for each plate we find the faintest targeted object in r-band magnitude with S/N larger than our S/N cut, with apparent magnitude r_cut (if this object is fainter than the nominal limit r_max for bright or faint plates, we set r_cut equal to this limit; r_max = 17.8 mag for bright plates and 20.2 mag for faint plates), and then assume that the selection function for that plate is given by a hyperbolic-tangent cutoff, centered on r_cut − 0.1 mag, and with a width parameter whose natural logarithm is −3 (≈0.05 mag), such that the total width of the cutoff is about 0.2 mag and the faintest object on the plate is about 2 widths from the center of the cutoff. The function value at the bright end is equal to the number of spectroscopic objects brighter than r_cut divided by the total number of targets brighter than r_cut. Thus, the plate-dependent selection function is given by

$$\begin{eqnarray} &&S(\mathrm{plate},r,g-r)\nonumber\\ &&\quad= \frac{\rm{No.\ of }\ \mathrm{spectroscopic\ objects\ with}\ {\it r}_{\mathrm{min}} \le {\it r} \le {\it r}_{\mathrm{cut}}}{\rm{No.\ of\ } \mathrm{targets\ with}\ {\it r}_{\mathrm{min}} \le {\it r} \le {\it r}_{\mathrm{cut}}} \nonumber\\ &&\qquad\times \left[1-\tanh \left(\frac{{\it r}-{\it r}_{\mathrm{cut}}+0.1}{\exp \left(-3\right)}\right)\right]\Big / 2\,, \end{eqnarray} \tag{ A2 }$$

where the numbers of objects are evaluated within the ≈7 deg² area of the plate in question and in the 0.48 ⩽ g − r ⩽ 0.55 G-star color range; r_min is 14.5 mag for bright plates and 17.8 mag for faint plates. The selection function is zero outside of the apparent r-band magnitude range of the plate ([14.5, 17.8] for bright plates and [17.8, 20.2] for faint plates).

We use this model both for the bright plates and the faint plates, although most bright plates are in fact consistent with being complete up to 17.8 mag. Figure 10 shows the distribution of Kolmogorov–Smirnov (K-S) probabilities that the spectroscopic sample for any given plate was selected from the target sample with this model for the selection function. All but seven plates have probabilities larger than 0.001 and the distribution of probabilities is relatively flat, as expected.

Zoom In Zoom Out Reset image size

Rather than using a smooth hyperbolic-tangent cutoff, we also tried a sharp cut at r_cut. With this model for the selection function, 79 plates have a K-S probability <0.05 (≈25% of the number of plates), as opposed to 30 plates in the hyperbolic-tangent-cutoff model (≈9% of the sample). Therefore, the smooth cutoff is necessary to fully model the selection function. The fact that the distribution of K-S probabilities in Figure 10 is not entirely flat is due to remaining details in the faint cutoff of the selection function, as we know that the selection function is flat at brighter magnitudes. This does not impact our analysis greatly, as most stars are much brighter than the cutoff (as compared to the scale over which the selection function changes near the cutoff).

The selection function is simplest in its native coordinates, survey plate, and r-band magnitude. For each value of g − r and [Fe/H], the r-dependent selection function above translates into a (different) spatial selection function through the use of the photometric-distance relation. The selection function projected into spatial coordinates for a typical value of g − r and [Fe/H] is shown in Figure 11. Near |b| = 90° the spectroscopic sample is relatively complete, whereas near the Galactic plane the selection is much less complete.

Zoom In Zoom Out Reset image size

We have posted a Python code that implements this model for the SEGUE selection function. It is publicly available at https://github.com/jobovy/segueSelect .

Figure 14 compares the observed distribution of vertical heights |Z| of the α-old G-dwarf sample to that predicted by the best-fit model. This prediction is obtained by running the best-fit density model integrated over the color–metallicity distribution through our model for the SEGUE selection function. There are 35 stars with magnitudes that should put them on faint plates, but that were observed on bright plates, and 551 stars in the opposite situation are cut from the data sample to show this comparison. The comparison between the data and the model is shown for all plates and for the bright and faint plates separately. The agreement between the model and observed distribution is excellent for all of these. Figure 15 shows a similar comparison for the distribution of Galactocentric radii of the data and in the model. The model correctly predicts the observed star counts for most Galactocentric radii, except the smallest around 5 kpc, where the model slightly overpredicts the number of stars (here, and in further comparisons below, the model around 8 kpc behaves somewhat erratically, as this is the boundary between 90° ⩽ l ⩽ 270° and −90° < l < 90° plates, and we do not use the finite extent of the plate in our model distributions). Also shown in this figure and in Figure 14 are models that only differ from the best-fit model in their radial scale length: a model with a scale length of 3 kpc and one with a scale length of 4 kpc. It is clear that these longer scale lengths are strongly ruled out by the data, as they strongly overpredict the star counts at large Galactocentric radii.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 1 lists the best-fit parameters for fits that only use (1) bright or faint plates, (2) b > 0° or b < 0° plates, or (3) |b| > 45° or |b| < 45° plates. The results from all of these different samples are roughly consistent with each other; we note that we can even measure the radial scale length with high-latitude plates (|b| > 45°) alone.

Figures 14 and 15 show comparisons between the observed star counts and the model, when we split the α-old sample into more metal-poor and more metal-rich sub-samples by cutting the sample at [Fe/H] = −0.7. Comparisons for when we split the α-old sample into two bins in [α/Fe], by cutting the sample at [α/Fe] = 0.35, are shown in Figure 16.

Zoom In Zoom Out Reset image size

Figure 17 compares the best-fit model to the observed star counts as a function of vertical height and Figure 18 shows this comparison as a function of Galactocentric radius, again removing 4 stars with magnitudes that should put them on faint plates but that were observed on bright plates and removing 43 stars in the opposite situation. We also show models whose parameters are the same as those of the best-fit model, but with shorter scale lengths of 2 and 3 kpc. The faint plates, which only contain 6% of the α-young sample, rule out a short scale length of 2 kpc for the α-young disk. The best-fit model provides a good fit to the observed star counts.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We again also list the best-fit parameters for fits that only use (1) bright or faint plates, (2) b > 0° or b < 0° plates, or (3) |b| > 45° or |b| < 45° plates in Table 2. The results from all of these different samples are again roughly consistent, except for the faint plates fit, which prefers even longer radial scale lengths, but faint plates only contain 6% of the α-young sample. We have also run fits for the α-young sample where we (1) employ a more conservative S/N cut of S/N > 30, (2) enlarge our sample with a less conservative S/N cut of S/N > 10, (3) remove stars on plates whose K-S probability for the spectroscopic sample to have been drawn from the underlying photometric sample combined with our model for the SEGUE selection function (see Figure 10) is smaller than 0.1, (4) use stars from the SEGUE database that were explicitly targeted as G-type stars (with all other log g, S/N, E(B − V) cuts), (5) remove stars with magnitudes that should put them on SEGUE bright plates, but that were observed as part of a faint plate and vice versa, and (6) artificially shift the metallicity distribution 0.1 dex toward the more metal-rich end. The results from these fits are all consistent with those obtained for our nominal sample with fiducial cuts.

Figures 17 and 18 show comparisons between the observed and predicted star counts for the α-young samples that are more metal-poor than the nominal α-young sample. The fit for the −0.6 < [Fe/H] <−0.3 sample is good, while the fit for the most metal-poor α-young sub-sample is not entirely satisfactory. Comparisons between the observed star counts and the model when we split the α-young sample into two, by cutting at [α/Fe] = 0.15, are shown in the lower two rows of Figure 16.