The Chemical Abundance Structure of the Inner Milky Way: A Signature of "Upside-down" Disk Formation

Stars in the vicinity of the Sun display a double-exponential distribution of heights $| z|$ above the disk midplane, often described as the superposition of a "thin disk" and "thick disk" (Gilmore & Reid 1983; Jurić et al. 2008). Young stars have cool kinematics that confine them to the thin disk, while stars with hotter kinematics are, on average, older, less metal-rich, and enhanced in the α-elements produced by core-collapse supernovae. Abundance studies show that the distribution of [α/Fe] is bimodal at solar or sub-solar [Fe/H], allowing chemical definition of "thin" and "thick" disk populations that display distinct kinematics (Bensby et al. 2003; Lee et al. 2011). A variety of mechanisms have been proposed for the origin of the thick disk, including debris of an accreted satellite (Abadi et al. 2003), heating of a pre-existing disk by satellites (Villalobos & Helmi 2008), gas-rich mergers or a turbulent gas-rich phase in the early star-forming disk (Brook et al. 2004; Bournaud et al. 2009; Forbes et al. 2012), nested flares of mono-age stellar populations (Minchev et al. 2015), or radial migration of stars outward from the kinematically hot inner disk (Schönrich & Binney 2009b). Two basic questions about the disk are: (1) were thick disk stars formed in a thick layer or heated over time from a thin distribution, and (2) was the formation of the thick disk a discrete event or extended over many orbital times?

In support of the "continuous" view, Bovy et al. (2012) show that mono-abundance populations in the disk, each defined by a narrow bin of $[\mathrm{Fe}/{\rm{H}}]$ and $[\alpha /\mathrm{Fe}]$, individually have exponential radial and vertical profiles. The superposition of multiple populations with different scale-heights gives rise to the double-exponential vertical profile observed for the full disk population and accounts for the correlation between [α/Fe] and kinematics. Inspired by this finding, Bird et al. (2013) show that similar behavior arises for mono-age populations in a high-resolution hydrodynamic simulation of the formation of a Milky Way-like galaxy, starting from cosmological initial conditions. Furthermore, Bird et al. (2013) show that the oldest disk stars are born in a vertically thick distribution from a turbulently supported interstellar medium (ISM), in accordance with previous numerical and analytic arguments about the influence of gravitational instability and feedback on the scale-height of the star-forming gas disk (Bournaud et al. 2009; Forbes et al. 2012). They dub this scenario "upside-down" disk formation, by analogy to "inside-out" growth in which the disk scale length increases with time (Kepner 1999; Pilkington et al. 2012). While subsequent heating does increase the scale-height of older populations (Bird et al. 2013, see their Figure 19, and J. Bird et al. 2017, in preparation), cosmological simulations suggest that the distinction between the thick and thin disks was largely imprinted at birth.

In this paper, we show that the combination of upside-down disk formation with a simple chemical evolution model can naturally explain the observed [α/Fe]–[Fe/H] distributions and metallicity distribution functions (MDFs) of the inner disk. The APOGEE survey (Majewski et al. 2017) has allowed mapping of the [α/Fe]–[Fe/H] distribution over much of the Milky Way disk through high-resolution near-infrared spectroscopy of more than 100,000 evolved stars. In most regions of the Galaxy, as in the solar neighborhood, stars lie on two distinct sequences in [α/Fe]–[Fe/H], though the locus of high-α stars appears to be independent of position within the disk (Hayden et al. 2015, hereafter H15). However, in the innermost radial zone mapped by H15, Galactocentric radius $3\,\mathrm{kpc}\lt R\lt 5\,\mathrm{kpc}$, stars lie along a sequence that resembles the evolutionary track of a simple one-zone chemical evolution model (e.g., Talbot & Arnett 1971; Tinsley 1980). In such a model, low metallicity stars occupy a plateau in [α/Fe] that reflects the characteristic yields of core-collapse supernovae, averaged over the stellar initial mass function (IMF). Once Type Ia supernovae become an important source of enrichment, they further increase [Fe/H], but they drive down [α/Fe] because they produce mainly iron-peak elements. The location of the "knee", defined here as the point at which the track has fallen by 0.1 dex in [α/Fe] from the initial plateau, depends mainly on the star formation efficiency (SFE). The value of [Fe/H] at late times depends mainly on the efficiency of galactic outflows, which controls the equilibrium between continuing enrichment from supernovae and loss of metals to stars and winds (Andrews et al. 2017; Weinberg et al. 2017, hereafter WAF). At $R=3\mbox{--}5\,\mathrm{kpc}$, most stars far from the midplane have enhanced [α/Fe] and sub-solar [Fe/H], while most stars near the midplane have $[\alpha /\mathrm{Fe}]\approx 0$ and super-solar [Fe/H]. We will show that this distribution can be readily explained by one-zone chemical evolution in a vertically contracting disk.

Explaining the bimodal [α/Fe]–[Fe/H] distribution seen at larger R requires something more complex than the chemical evolution models considered here, such as mixing of stellar populations by radial migration (Sellwood & Binney 2002; Schönrich & Binney 2009a), sharp transitions in the Galaxy accretion and star formation history (Chiappini et al. 1997), or an increase of outflow efficiency at late times that induces reverse evolution toward lower [Fe/H] (WAF). Here, we focus on the inner disk, where simple models appear applicable. We suspect that the more complex structure of the outer disk arises because radial migration has a larger impact there, while the inner disk is dominated by stars that remain close to their birth radius. Tentative support for this view comes from the modeling of MDFs presented by H15 (see their Figure 10). A full model for the formation and evolution of the Milky Way disk must account for its full spatial and chemo-dynamical structure and explain the distributions of elements that trace different nucleosynthetic pathways. This paper presents idealized models intended to highlight one aspect of this evolution.

Our data sample is a subset of that studied by H15, drawn from the APOGEE survey (Majewski et al. 2017) of the Sloan Digital Sky Survey III (SDSS-III, Eisenstein et al. 2011). It consists of red giant and sub-giant stars with elemental abundances inferred from infrared (H-band) spectroscopy with spectral resolution $\lambda /{\rm{\Delta }}\lambda \approx \mathrm{20,000}$ (Nidever et al. 2015). Elemental abundances are inferred from the APOGEE Stellar Parameters and Chemical Abundances Pipeline (ASPCAP; García Pérez et al. 2016). In brief, the parameters of each star are determined by ${\chi }^{2}$ minimization across the full spectrum, using 6-parameter model grids: T_eff, $\mathrm{log}g$, an overall metallicity [M/H], scalings of α-elements relative to the overall metallicity [α/M], and individually varied abundances of C and N. With free variations of α, C, and N, the metallicity M is responsive mainly to iron-peak elements. We (like H15) refer to [α/M] and [M/H] by the more conventional labels [α/Fe] and [Fe/H]. We apply the same cuts as H15, summarized here in Table 1, to select main sample, evolved stars with high-quality abundance determinations. These abundances are the ones provided in SDSS DR12 (Alam et al. 2015; Holtzman et al. 2015). Further details of APOGEE field and target selection can be found in Zasowski et al. (2013).

Table 1.  Cuts Made on APOGEE Data for Sample Selection

Cuts	Notes
$3\,\mathrm{kpc}\lt R\lt 5\,\mathrm{kpc}$	Inner annulus of H15
$1.0\lt \mathrm{log}g\lt 3.8$	Choose giants only
$3500\,{\rm{K}}\lt {T}_{\mathrm{eff}}\lt 5500\,{\rm{K}}$	High-quality abundances
${\rm{S}}/{\rm{N}}\gt 80$	"
APOGEE_TARGET1 bits ∈ 11,12,13	Main survey targets only
ASPCAPFLAG bits ∉ 23	Remove bad stars

Download table as:  ASCIITypeset image

H15 estimate distances from spectroscopically determined stellar parameters, PARSEC isochrones (Bressan et al. 2012), and extinction-corrected apparent magnitudes (see H15 and Holtzman et al. 2015 for details). They use these distances to assign stars to bins in Galactocentric radius R and vertical distance from the midplane $| z|$. In this paper, we concentrate on H15's innermost radial bin 3 kpc $\lt \,R\,\lt$ 5 kpc, where stars lie along what is plausibly a single evolutionary track. After the cuts listed in Table 1, we are left with 3753 stars: 2435 at $| z| =0\mbox{--}0.5\,\mathrm{kpc}$, 849 at $| z| =0.5\mbox{--}1.0\,\mathrm{kpc}$, and 469 at $| z| =1.0\mbox{--}2.0\,\mathrm{kpc}$. The top panels of Figure 1 plot the MDFs of stars in these three vertical layers, and the bottom panels plot the distribution of stars in the [α/Fe]–[Fe/H] plane, together with the evolutionary track described in the following section. Equivalent information in different formats appears in Figures 4 and 5 of H15. In the main body of the paper, we treat these measured MDFs as representing the fraction of stars born at each metallicity, even though APOGEE observes evolved stars in particular magnitude ranges. In Appendix B, we show that the impact of APOGEE selection on the MDF is small, at least given the age–metallicity relation implied by our models, justifying this approximate treatment.

Zoom In Zoom Out Reset image size

We adopt oxygen as the representative α-element in our theoretical model. The H15 measurements depend on several α-elements with oxygen dominant among them. Within our model, the defining characteristic of α elements is that they are produced by core-collapse supernovae only, with a yield that is independent of metallicity.

We calculate chemical evolution within a single galactic zone according to the following specifications.

• 1.  
  The galactic zone consists of a supply of gas out of which stars are formed over time. This gas is initially composed of pristine hydrogen and helium, and it is homogeneously mixed at all times. Any new metals introduced by supernovae are immediately mixed into the existing gas supply.
• 2.  
  Stars are formed from the gas supply with a constant efficiency ${\dot{M}}_{* }/{M}_{g}$, which we specify in terms of a SFE timescale ${\tau }_{* }={M}_{g}/{\dot{M}}_{* }$. We choose ${\tau }_{* }=2.0\,\mathrm{Gyr}$, based on the observationally inferred timescale for molecular gas in local disk galaxies (Leroy et al. 2008).
• 3.  
  The star formation rate (SFR) evolves in a linear-exponential fashion, i.e.,

  $$\begin{eqnarray}&&{\dot{M}}_{* }(t)={\dot{M}}_{* }({t}_{\mathrm{sf}})\times \left(\displaystyle \frac{t}{{t}_{\mathrm{sf}}}\right){e}^{-t/{t}_{\mathrm{sf}}}\,,\end{eqnarray} \tag{ 1 }$$

  where ${t}_{\mathrm{sf}}$ is the star formation timescale. Together with assumption 2, this implies a linear-exponential form for the gas supply ${M}_{g}(t)={\tau }_{* }{\dot{M}}_{* }(t)$. We consider several different values for ${t}_{\mathrm{sf}}$. Appendix A shows results for an exponential star formation history, but the linear-exponential model corresponds better to the predictions of galaxy formation simulations (Simha et al. 2014), and it yields a somewhat better fit to the H15 data.
• 4.  
  Outflow of gas from the zone, due to stellar feedback, occurs at a rate ${\dot{M}}_{\mathrm{out}}$. We parameterize outflow efficiency with the mass-loading factor η, such that $\eta ={\dot{M}}_{\mathrm{out}}/{\dot{M}}_{* }$. We adjust the value of η depending on ${t}_{\mathrm{sf}}$ to produce similar peaks in the [Fe/H] distribution (see discussion below).
• 5.  
  Each population of stars forms with masses distributed according to a Kroupa (2001) IMF, which prescribes a double power-law form for

  $$\begin{eqnarray}&&\zeta (M)\propto \displaystyle \frac{{dN}(M)}{{dM}}\propto {M}^{-\alpha }\end{eqnarray} \tag{ 2 }$$

  where

  $$\begin{eqnarray}\alpha =\left\{\begin{array}{ll}1.8, & 0.08\lt M/{M}_{\odot }\lt 0.5\\ 2.7, & 0.5\lt M/{M}_{\odot }\lt 1.0\\ 2.3, & 1.0\lt M/{M}_{\odot }\end{array}\right..\end{eqnarray} \tag{ 3 }$$

  We truncate the IMF at ${M}_{\min }=0.08\,{M}_{\odot }$ and ${M}_{\max }\,=100\,{M}_{\odot }$.
• 6.  
  There are two sources of new metal production: core-collapse supernovae (CCSNe) and Type Ia supernovae (SN Ia). We define the yield parameters ${m}_{{\rm{O}}}^{\mathrm{CC}}$, the mass of oxygen produced by CCSNe per unit mass of star formation, ${m}_{\mathrm{Fe}}^{\mathrm{CC}}$, the corresponding quantity for iron, and ${m}_{\mathrm{Fe}}^{\mathrm{Ia}}$, the mass of iron produced by SN Ia over a time interval 12.5 Gyr per unit mass of star formation. We adopt ${m}_{{\rm{O}}}^{\mathrm{CC}}=0.015$ and ${m}_{\mathrm{Fe}}^{\mathrm{CC}}=0.0012$, based on a Kroupa IMF and the CCSN yields of Chieffi & Limongi (2004) and Limongi & Chieffi (2006). We adopt ${m}_{\mathrm{Fe}}^{\mathrm{Ia}}=0.0017$ based on the SN Ia rate of Maoz et al. (2012) and an iron yield of $0.77\,{M}_{\odot }$ (Iwamoto et al. 1999) per supernova (see Andrews et al. 2017 for more detailed discussion). These yields are defined to be net quantities: they do not include the amount of oxygen and iron incorporated into the stellar population from the ISM at birth, but only the newly produced amounts of these elements.
• 7.  
  CCSNe go off instantaneously upon the birth of a stellar population; this is a reasonable approximation because we adopt star formation histories that are smooth on much longer timescales. However, it does assume that the "lock up" of metals in the warm or hot phases of the ISM (Schönrich & Binney 2009a) can be ignored.
• 8.  
  SN Ia follow a power-law delay time distribution (DTD) after a minimum delay time t_d:

  $$\begin{eqnarray}R(t)=\left\{\begin{array}{ll}{R}_{0}{t}^{-1.1} & t\gt {t}_{d}\\ 0 & t\lt {t}_{d}.\end{array}\right.\end{eqnarray} \tag{ 4 }$$

  The form is based on Maoz et al. (2012), and the normalization R₀ cancels out of our calculations once we have specified the iron yield ${m}_{\mathrm{Fe}}^{\mathrm{Ia}}$. We adopt ${t}_{d}=0.15\,\mathrm{Gyr}$. We have tested that the choice of minimum delay time has little impact on our results, as the analytic framework of WAF would suggest (see, e.g., their Figures 4 and 11).
• 9.  
  A star of mass M lives for a time ${t}_{* }=10\,\mathrm{Gyr}\,\times {(M/{M}_{\odot })}^{-3.5}$, then returns to the ISM as a mass of gas M − ${M}_{\mathrm{rem}}$ with the same chemical composition the star had at birth. (Stars above $8\,{M}_{\odot }$ additionally return elements with the CCSN net yields described above.) The remnant mass as a function of initial mass is

  $$\begin{eqnarray}{M}_{\mathrm{rem}}(M)=\left\{\begin{array}{ll}0.394\,{M}_{\odot }+0.109M & M\lt 8\,{M}_{\odot }\\ 1.44\,{M}_{\odot } & M\gt 8\,{M}_{\odot }.\end{array}\right.\end{eqnarray} \tag{ 5 }$$

  (see Kalirai et al. 2008). The lifetime formula can be inverted to yield the main sequence turnoff mass as a function of time:

  $$\begin{eqnarray}&&{M}_{\mathrm{to}}(t)={\left(\displaystyle \frac{t}{10\mathrm{Gyr}}\right)}^{-1/3.5}{M}_{\odot }.\end{eqnarray} \tag{ 6 }$$

With the assumptions of point 9, the fraction of a mass of stars formed in a short interval at time $t^{\prime}$ that has been recycled to the ISM by time t is

$$\begin{eqnarray}&&{F}_{\mathrm{rec}}(t,t^{\prime} )=\displaystyle \frac{{\int }_{{M}_{\mathrm{to}}(t-t^{\prime} )}^{{M}_{\max }}M\zeta (M)[1-{f}_{\mathrm{rem}}(M)]\,{dM}}{{\int }_{{M}_{\min }}^{{M}_{\max }}M\zeta (M){dM}}\,,\end{eqnarray} \tag{ 7 }$$

where ${f}_{\mathrm{rem}}(M)={M}_{\mathrm{rem}}(M)/M$ is the fraction of a star of initial mass M that is left in the stellar remnant. While our lifetime formula is inaccurate at high masses, it adequately captures the important behavior, which is that high mass stars recycle their envelopes quickly compared to the timescales for SN Ia enrichment and variation of the SFR. For a Kroupa IMF, the recycled mass fraction is $r={F}_{\mathrm{rec}}(t,t^{\prime} )=0.37$, 0.40, and 0.45 for $t-t^{\prime} =1$, 2, and 10 Gyr, respectively. The difference in chemical evolution between our complete time-dependent recycling and an instantaneous approximation with r = 0.4 is small (see Figure 7 of WAF).

With these assumptions in place, we may write down the differential equations governing the amount of oxygen and iron in the ISM. Metals are injected into the ISM by CCSNe, SN Ia, and recycling, and they are removed from the ISM by star formation and outflows. For oxygen,

$$\begin{eqnarray}{\dot{M}}_{{\rm{O}}} & = & {m}_{{\rm{O}}}^{\mathrm{CC}}{\dot{M}}_{* }-{\dot{M}}_{* }{Z}_{{\rm{O}}}-\eta {\dot{M}}_{* }{Z}_{{\rm{O}}}\\ & & +{\displaystyle \int }_{0}^{t}{\dot{M}}_{* }(t^{\prime} ){Z}_{{\rm{O}}}(t^{\prime} ){\dot{F}}_{\mathrm{rec}}(t,t^{\prime} ){dt}^{\prime} \\ & = & {m}_{{\rm{O}}}^{\mathrm{CC}}\displaystyle \frac{{M}_{g}}{{\tau }_{* }}-(1+\eta )\displaystyle \frac{{M}_{{\rm{O}}}}{{\tau }_{* }}\\ & & +{\displaystyle \int }_{0}^{t}\displaystyle \frac{{M}_{{\rm{O}}}(t^{\prime} )}{{\tau }_{* }}{\dot{F}}_{\mathrm{rec}}(t,t^{\prime} ){dt}^{\prime} ,\end{eqnarray} \tag{ 8 }$$

where ${Z}_{{\rm{O}}}={M}_{{\rm{O}}}/{M}_{g}$ is the oxygen mass fraction and ${\dot{F}}_{\mathrm{rec}}(t,t^{\prime} )$ is the derivative of Equation (7) with respect to t. For iron, we have the additional contribution of SN Ia:

$$\begin{eqnarray}{\dot{M}}_{\mathrm{Fe}} & = & {m}_{\mathrm{Fe}}^{\mathrm{CC}}\displaystyle \frac{{M}_{g}}{{\tau }_{* }}-(1+\eta )\displaystyle \frac{{M}_{\mathrm{Fe}}}{{\tau }_{* }}+{\displaystyle \int }_{0}^{t}\displaystyle \frac{{M}_{\mathrm{Fe}}(t^{\prime} )}{{\tau }_{* }}{\dot{F}}_{\mathrm{rec}}(t,t^{\prime} ){dt}^{\prime} \\ & & +{m}_{\mathrm{Fe}}^{\mathrm{Ia}}\displaystyle \frac{{\displaystyle \int }_{0}^{t}R(t-t^{\prime} ){M}_{g}(t^{\prime} )/{\tau }_{* }\,{dt}^{\prime} }{{\displaystyle \int }_{0}^{\infty }R(t^{\prime} ){dt}^{\prime} }.\end{eqnarray} \tag{ 9 }$$

We solve these equations for oxygen and iron evolution with a numerical code that uses simple Euler integration to evolve the zone from t = 0 to $t=12.5\,\mathrm{Gyr}$, with a step size of ${\rm{\Delta }}t=10\,\mathrm{Myr};$ that is,

$$\begin{eqnarray}&&{M}_{{\rm{X}},i}={M}_{{\rm{X}},i-1}+{\dot{M}}_{{\rm{X}},i-1}{\rm{\Delta }}t.\end{eqnarray} \tag{ 10 }$$

We evaluate the recycling term from the cumulative values of ${F}_{\mathrm{rec}}(t,t^{\prime} )$ at the beginning and end of the step. We prescribe a star formation history rather than an infall history (see Vincenzo et al. (2017) for a mathematical discussion of the relation between these), and the gas mass is therefore implicitly specified by ${M}_{g}(t)={\tau }_{* }{\dot{M}}_{* }(t)$, with no need to integrate it separately.

Our assumptions are similar to those of many one-zone chemical evolution models, and they closely follow those adopted for analytic modeling by WAF. WAF show that when the integral recycling terms in Equations (8) and (9) are replaced by the instantaneous approximation ${{rM}}_{X}(t)/{\tau }_{* }$, then the resulting equations for oxygen and iron evolution can be solved analytically if one adopts an exponential DTD for SN Ia enrichment and a constant, exponential, or linear-exponential star formation history. These solutions show that for a wide range of model parameters the values of $[\mathrm{Fe}/{\rm{H}}]$ and $[{\rm{O}}/\mathrm{Fe}]$ converge within a few Gyr to "equilibrium" values at which sources and sinks balance, reproducing numerical results (Andrews et al. 2017). The equilibrium $[\mathrm{Fe}/{\rm{H}}]$ depends mainly on yields and on the value of η, while the equilibrium $[{\rm{O}}/\mathrm{Fe}]$ depends mainly on yields. The importance of equilibrium, and the key role of outflow efficiency in regulating abundances, has been highlighted in theoretical studies of the galaxy mass–metallicity relation (Finlator & Davé 2008; Peeples & Shankar 2011; Davé et al. 2012).

In the cosmological simulation of Bird et al. (2013), the median height of the star-forming gas layer decreases from about 0.8 kpc to about 0.2 kpc between $t\approx 2\,\mathrm{Gyr}$ and $t\approx 6\,\mathrm{Gyr}$, as the ratio of rotational velocity to turbulent vertical velocity grows rapidly. After that, the gas layer continues to get thinner, but at a much slower rate. Rather than mimic a single simulation, we adopt an idealized model in which the scale-height z_h of newly formed stars decreases linearly from 0.8 to 0.2 kpc over a collapse timescale t_c and remains constant thereafter:

$$\begin{eqnarray}{z}_{h}(t,{t}_{c})=\left\{\begin{array}{ll}0.8\,\mathrm{kpc}-0.6\,\mathrm{kpc}\times (t/{t}_{c}) & t\leqslant {t}_{c},\\ 0.2\,\mathrm{kpc} & t\gt {t}_{c}.\end{array}\right.\end{eqnarray} \tag{ 11 }$$

Figure 2 compares these scale-height histories to the star formation histories for the values of t_c and ${t}_{\mathrm{sf}}$ that we consider below. In a model where gravitationally driven turbulence of a gas-rich disk governs both the thickness of the gas layer and the efficiency of star formation (Bournaud et al. 2009; Forbes et al. 2012), we might expect the collapse time to be linked to the evolution of the SFR. We do not impose this condition a priori, but we find that our fits to the H15 data prefer comparable collapse and star formation timescales.

Zoom In Zoom Out Reset image size

We evolve our model for 12.5 Gyr. As noted previously, we choose the gas supply so that the SFR is

$$\begin{eqnarray}&&\displaystyle \frac{{{dN}}_{* }}{{dt}}\propto \displaystyle \frac{{M}_{g}(t)}{{\tau }_{* }}\propto {{te}}^{-t/{t}_{\mathrm{sf}}}.\end{eqnarray} \tag{ 12 }$$

The ${\rm{\Delta }}{N}_{* }={\rm{\Delta }}t\times ({{dN}}_{* }/{dt})$ stars formed in a given timestep are assigned z values drawn from the distribution

$$\begin{eqnarray}&&h(z,{t}_{c})\propto {e}^{-z/{z}_{h}(t,{t}_{c})}\,,\end{eqnarray} \tag{ 13 }$$

where ${z}_{h}(t,{t}_{c})$ is given by Equation (11). Thus, for a given t_c, the fraction of stars formed at time t with ${z}_{1}\lt z\lt {z}_{2}$ is given by

$$\begin{eqnarray}&&\displaystyle \frac{{N}_{{z}_{1}\to {z}_{2}}}{{N}_{\mathrm{tot}}}=\displaystyle \frac{{\int }_{{z}_{1}}^{{z}_{2}}{e}^{-z/{z}_{h}(t,{t}_{c})}{dz}}{{\int }_{0}^{\infty }{e}^{-z/{z}_{h}(t,{t}_{c})}{dz}}.\end{eqnarray} \tag{ 14 }$$

We ignore any subsequent changes to the vertical structure of the disk, e.g., from dynamical heating or adiabatic compression, a point we return to at the end of Section 4. Each star formed is assigned oxygen and iron abundances according to the abundances in the ISM at time t, which are determined by the chemical evolution model described above. Scatter is added to these values, drawn from a Gaussian distribution with $\sigma =0.05$ dex. We implement scatter partly to enable better visual interpretation of our plots and partly to mimic effects of observational errors, although our choice of value for the characteristic uncertainty somewhat exaggerates the scatter at high metallicities.

Figure 3 displays the predicted distribution of stars in [O/Fe]–$[\mathrm{Fe}/{\rm{H}}]$, in the form of color-coded two-dimensional (2D) histograms, in the same three vertical layers used in Figure 1. We consider t_c values of 1.0, 2.5, and 4.5 Gyr and ${t}_{\mathrm{sf}}$ values of 1.0, 3.0, and 6.0 Gyr (see Figure 2), making nine models in total. In each panel, the central track shows the model abundances calculated without scatter, and colored points mark the abundance values at t = 0.1, 0.5, 1.0, 2.0, 5.0, and 10.0 Gyr. The tracks and colored points are identical in a given column of the plot because z and t_c do not enter into the chemical evolution. As explained by WAF, equilibrium abundances depend mainly on yields and η if ${t}_{\mathrm{sf}}$ is long, but as ${t}_{\mathrm{sf}}$ approaches the gas depletion time ${\tau }_{* }/(1+\eta -r)$ the equilibrium abundance grows, and the approach to equilibrium gets slower. For ${t}_{\mathrm{sf}}=6.0\,\mathrm{Gyr}$ and $3.0\,\mathrm{Gyr}$, we have adopted $\eta =1.17$ and 1.50, so that both models yield $[\mathrm{Fe}/{\rm{H}}]\approx +0.25$ at $t=10\,\mathrm{Gyr}$. For ${t}_{\mathrm{sf}}=1.0\,\mathrm{Gyr}$, we choose $\eta =2.83$, which gives $[\mathrm{Fe}/{\rm{H}}]\approx +0.25$ at $t=5\,\mathrm{Gyr}$, but because of the rapid truncation of the SFR, most stars form below this metallicity. Tracks and time-stamped abundances are fairly similar in the ${t}_{\mathrm{sf}}=6.0$ and $3.0\,\mathrm{Gyr}$ models, though one can see the slightly slower approach to equilibrium for ${t}_{\mathrm{sf}}=3.0\,\mathrm{Gyr}$. The evolutionary track and, especially, abundance versus time are significantly different for ${t}_{\mathrm{sf}}=1.0\,\mathrm{Gyr}$.

Zoom In Zoom Out Reset image size

As shown previously in Figure 1, the evolutionary track of our ${t}_{\mathrm{sf}}=3.0\,\mathrm{Gyr}$ model follows the observed locus of H15 stars quite well. This agreement indicates that, given our adopted CCSN and SN Ia yields, our adopted values of ${\tau }_{* }$ and η are reasonable. We note in particular that a much higher (lower) SFE would drive the knee of the evolutionary track to higher (lower) $[\mathrm{Fe}/{\rm{H}}]$, in disagreement with the data, though we have not investigated the degree to which a change in ${\tau }_{* }$ could be compensated by changes in the star formation history or SN Ia DTD. The predicted $[{\rm{O}}/\mathrm{Fe}]$ at high $[\mathrm{Fe}/{\rm{H}}]$ is slightly sub-solar, while the APOGEE measurements yield solar $[\alpha /\mathrm{Fe}]$. Given the uncertainties in observational calibrations, solar abundance values, and supernova yields, we do not ascribe much significance to this discrepancy, but models with longer ${t}_{\mathrm{sf}}$ and thus higher equilibrium $[{\rm{O}}/\mathrm{Fe}]$ would fit the data better with other assumptions held fixed. Physically, a more rapidly declining star formation history leads to lower $[{\rm{O}}/\mathrm{Fe}]$ at late times because the ratio of CCSNe to SN Ia is lower.

The choices of t_c and ${t}_{\mathrm{sf}}$ affect the distribution of stars, both along the tracks and among layers of the disk. Examining the middle column isolates the effect of t_c. The most distinctive impact can be seen in the $| z| =1\mbox{--}2\,\mathrm{kpc}$ layers of the t_c = 1.0 and $2.5\,\mathrm{Gyr}$ models. For the shorter collapse timescale, the distribution of stars in this layer is bimodal, with one group of stars near the knee of the track (formation times centered at $t\approx 0.5\,\mathrm{Gyr}$) and a second group near the end of the track (formation times centered at $t\approx 5\,\mathrm{Gyr}$). The first group forms during disk collapse, when the scale-height is still fairly high. By $t=1.0\,\mathrm{Gyr}$ the scale-height has dropped to 0.2 kpc, and although the SFR between t = 1 and 3 Gyr is high, only a modest fraction of all stars are formed in this interval (see Figure 2), and a small fraction of these have $| z| \gt 1\,\mathrm{kpc}$. The disk remains thin at later times, but more than 60% of the stars form between 3 and 10 Gyr, and the high-$| z|$ tail of this population is sufficient to produce the second group of stars, at high $[\mathrm{Fe}/{\rm{H}}]$. The central model, with ${t}_{c}=2.5\,\mathrm{Gyr}$, still has a bimodal distribution of stars in the high-$| z|$ later, but now the α-enhanced low metallicity group is far more prominent because the disk collapse is slower, and more stars form while the scale-height remains high.

Other differences in the middle column of Figure 3 are more subtle but follow similar trends. In the high-$| z|$ layer, the mean metallicity of the low-$[\mathrm{Fe}/{\rm{H}}]$ group increases as t_c increases, as the population can evolve toward higher metallicity before the disk contracts to low scale-height. For ${t}_{c}=1.0\,\mathrm{Gyr}$, the mid-$| z|$ layer is dominated by stars formed at late times when the disk is thin, but for ${t}_{c}=4.5\,\mathrm{Gyr}$, this layer is dominated by intermediate-age stars formed during collapse. The central, ${t}_{c}=2.5\,\mathrm{Gyr}$ model has the broadest distribution of age and $[\mathrm{Fe}/{\rm{H}}]$ in the mid-$| z|$ layer because these two groups join to make a single continuous distribution. The low-$| z|$ layers appear similar in all three models because a large fraction of stars form after the disk has reached its 0.2 kpc scale-height.

Each 2D histogram in Figure 3 is normalized relative to the maximum density in the same layer, similar to the analogous figure of H15 (their Figure 4). However, the relative number of stars in the higher $| z|$ layers depends strongly on t_c and to some degree on ${t}_{\mathrm{sf}};$ e.g., in the central column, the fraction of stars with $| z| =1\mbox{--}2\,\mathrm{kpc}$ increases from 1% to 5% between the t_c = 1.0 and $4.5\,\mathrm{Gyr}$ models. The relative numbers of stars in the H15 sample are strongly shaped by the survey strategy and selection function of APOGEE, so we cannot use them as a model test. A measurement of the $| z|$-distribution of stars with $R=3\mbox{--}5\,\mathrm{kpc}$ provides an additional diagnostic for distinguishing models, a point we return to in the discussion of Figure 8 below.

The impact of ${t}_{\mathrm{sf}}$ can be seen by looking across the middle (${t}_{c}=2.5\,\mathrm{Gyr}$) row of Figure 3. A shorter ${t}_{\mathrm{sf}}=1.0\,\mathrm{Gyr}$ (left column) leads to a broader distribution of $[\mathrm{Fe}/{\rm{H}}]$ in the high-$| z|$ layer, as rapid early star formation allows greater chemical enrichment before disk contraction. The low-$| z|$ layer is populated mainly by stars formed between t = 1 and $5\,\mathrm{Gyr}$, but despite the narrower age range, the metallicity range (roughly $[\mathrm{Fe}/{\rm{H}}]=-0.5$ to $+0.25$) is broader than that of the corresponding layer for ${t}_{\mathrm{sf}}=3.0\,\mathrm{Gyr}$. Moving to ${t}_{\mathrm{sf}}=6.0\,\mathrm{Gyr}$ (right column) leads to a much stronger concentration of low-$| z|$ stars near the equilibrium metallicity, both because the approach to equilibrium is faster for longer ${t}_{\mathrm{sf}}$ and because this model forms a larger fraction of its stars at late times (see Figure 2). At high $| z|$, this model shows a bimodal distribution similar to that of the (${t}_{\mathrm{sf}}=3.0\,\mathrm{Gyr},{t}_{c}=1.0\,\mathrm{Gyr}$) model discussed previously, again with a low-$[\mathrm{Fe}/{\rm{H}}]$ group formed during disk collapse and a high-$[\mathrm{Fe}/{\rm{H}}]$ group that is the high-$| z|$ tail of the much larger number of stars formed at late times with a 0.2 kpc scale-height. Most of the models in Figure 3 show some degree of this bimodality at high $| z|$, but the relative importance of the two groups depends on the specific values of ${t}_{\mathrm{sf}}$ and t_c. More generally, while there is some tradeoff between these two parameters, their impact is not strongly degenerate once all three layers are considered.

Figure 4 plots MDFs, the fraction of stars in bins of $[\mathrm{Fe}/{\rm{H}}]$, for the same nine models and three $| z|$-layers shown in Figure 3. Because all stars in a given model lie along a single evolutionary track, the one-dimensional MDF together with the track itself retains all of the information from Figure 3. The MDF allows easier comparison to the H15 measurements, shown by the orange curves in Figure 4. The ${t}_{\mathrm{sf}}=1.0\,\mathrm{Gyr}$ models are a consistently poor match to the data in the low-$| z|$ layer, predicting MDFs that are too broad and too symmetric. Changes to η, supernova yields, or ${\tau }_{* }$ could shift the centroid of the MDF, but the broad and symmetric shape is a generic consequence of the short star formation timescale (WAF; Andrews et al. 2017). The observed MDF at low $| z|$ instead shows the peaked, negatively skewed form expected when continuing infall allows a large fraction of stars to form near the equilibrium abundance, and the ${t}_{\mathrm{sf}}=3.0\,\mathrm{Gyr}$ and ${t}_{\mathrm{sf}}=6.0\,\mathrm{Gyr}$ models reproduce it much better. The high-$| z|$ MDF of these models is nearly always bimodal, but the relative amplitudes of the low-$[\mathrm{Fe}/{\rm{H}}]$ and high-$[\mathrm{Fe}/{\rm{H}}]$ peaks varies strongly from one $({t}_{\mathrm{sf}},{t}_{c})$ combination to another. The data do not show a clear bimodality in this layer, and the $({t}_{\mathrm{sf}}=3.0\,\mathrm{Gyr},{t}_{c}=2.5\,\mathrm{Gyr})$ model, for which the two peaks merge into a continuous broad distribution, provides the best match to the data. This model and the $({t}_{\mathrm{sf}}\,=3.0\,\mathrm{Gyr},{t}_{c}=4.5\,\mathrm{Gyr})$ model also provide the best match to the $| z| =0.5\mbox{--}1\,\mathrm{kpc}$ MDF, for which models with shorter t_c or longer ${t}_{\mathrm{sf}}$ predict too many high-$[\mathrm{Fe}/{\rm{H}}]$ stars. Our best models reproduce the $| z|$-dependence of the MDF better than the simulation-based model of Miranda et al. (2016; see their Figure 8), though the trend toward a broader and more symmetric MDF at high $| z|$ is also present in their simulations, and the innermost annulus they consider is 7–8 kpc rather than the 3–5 kpc range studied here.

Zoom In Zoom Out Reset image size

To quantify the level of agreement seen in Figure 4, we compute the squared difference of the MDF histograms

$$\begin{eqnarray}&&{{\rm{\Delta }}}^{2}=\displaystyle \sum _{i=0}^{N}{({f}_{* ,i}^{\mathrm{data}}-{f}_{* ,i}^{\mathrm{model}})}^{2}\,,\end{eqnarray} \tag{ 15 }$$

where the sum is over the N = 75 $[\mathrm{Fe}/{\rm{H}}]$ bins. Figure 5 plots these values for the nine models in each of the three $| z|$ layers. In the low-$| z|$ layer, the ${t}_{\mathrm{sf}}=1.0\,\mathrm{Gyr}$ model performs very poorly, while the ${t}_{\mathrm{sf}}=3.0\,\mathrm{Gyr}$ and $6.0\,\mathrm{Gyr}$ models reproduce the data about equally well, with almost no dependence on t_c. In the mid-$| z|$ layer, the ${t}_{\mathrm{sf}}=6.0\,\mathrm{Gyr}$ model fares poorly, especially for low t_c, and the ${t}_{\mathrm{sf}}=3.0\,\mathrm{Gyr}$ model performs best at all t_c. The high-$| z|$ layer strongly disfavors ${t}_{c}=1\,\mathrm{Gyr}$, for any ${t}_{\mathrm{sf}}$, and it again disfavors the ${t}_{\mathrm{sf}}=6.0\,\mathrm{Gyr}$ models at all t_c. Overall, the ${t}_{\mathrm{sf}}=3.0\,\mathrm{Gyr}$ models with t_c = 2.5 or $4.5\,\mathrm{Gyr}$ are clearly the most successful, with ${t}_{c}=2.5\,\mathrm{Gyr}$ slightly preferred. We have not attempted to formally evaluate parameter uncertainties or a global goodness-of-fit, in part because our models seem too idealized to warrant such an approach. In particular, a global ${\chi }^{2}$ or maximum likelihood statistic would assign the most weight to the low-$| z|$ MDF, which is based on the most stars and thus has the smallest statistical errors, but given that our models are likely incomplete, the level of agreement in the mid-$| z|$ and high-$| z|$ layers is a better diagnostic of star formation and disk collapse timescales than small improvements in reproducing the low-$| z|$ data. A formal fitting procedure would also require accounting for systematic uncertainties in the abundance data and for the influence of APOGEE's evolved star selection on the MDF.

Zoom In Zoom Out Reset image size

Figure 6 plots the first four moments of the MDFs—mean, variance, skewness, and kurtosis—for the nine models and the data in each $| z|$-layer. Our precise definitions of these moments are given in Table 2. In each panel, the horizontal black line marks the measured moment for the data and the gray band marks the 1σ error on the mean, determined from 20 bootstrap resamplings of the original data set. Because the MDFs are sometimes strongly non-Gaussian (e.g., bimodal), even four moments do not fully characterize their shapes. Overall, the moments comparison also shows the ${t}_{\mathrm{sf}}=3.0\,\mathrm{Gyr}$ models with ${t}_{c}=2.5$ or $4.5\,\mathrm{Gyr}$ to be the most successful. However, the ${t}_{\mathrm{sf}}=6.0\,\mathrm{Gyr}$ models give better agreement with the skewness of the low-$| z|$ MDF, and the rejection of ${t}_{c}=1.0\,\mathrm{Gyr}$ models is much more resounding for the ${{\rm{\Delta }}}^{2}$ statistic than it is for MDF moments.

Zoom In Zoom Out Reset image size

Table 2.  Statistical Quantities Used to Compare Model MDFs to the H15 Data; x_i Is the Value of [Fe/H] in the ith Bin, N is the Number of Bins, ${N}_{* ,i}$ Is the Number of Stars in the ith Bin, and ${f}_{* ,i}={N}_{* ,i}/{N}_{* ,\mathrm{tot}}$ is the Fraction of Stars in the ith bin

Name	Equation
${{\rm{\Delta }}}^{2}$	${\sum }_{i=0}^{N}{({f}_{* ,i}^{\mathrm{data}}-{f}_{* ,i}^{\mathrm{model}})}^{2}$
Mean	$\bar{x}=\tfrac{1}{N}{\sum }_{i=0}^{N}{x}_{i}{N}_{* i}$
Variance	${\sigma }^{2}=\tfrac{1}{N}{\sum }_{i=0}^{N}{({x}_{i}-\bar{x})}^{2}{N}_{* i}$
Skewness	$\tfrac{1}{N}{\sum }_{i=0}^{N}{[({x}_{i}-\bar{x})/\sigma ]}^{3}{N}_{* i}$
Kurtosis	$\tfrac{1}{N}{\sum }_{i=0}^{N}{[({x}_{i}-\bar{x})/\sigma ]}^{4}{N}_{* i}-3$

Download table as:  ASCIITypeset image

Figure 7 returns to a more qualitative comparison between our central model (${t}_{\mathrm{sf}}=3.0\,\mathrm{Gyr}$, ${t}_{c}=2.5\,\mathrm{Gyr})$ and the H15 data. In each layer, we select a number of stars similar to that in the H15 data, and we plot them (including the 0.05 dex of added scatter in $[{\rm{O}}/\mathrm{Fe}]$ and $[\mathrm{Fe}/{\rm{H}}]$) in the same manner as the lower three panels of Figure 1, which we repeat here to allow direct visual comparison. The model effectively reproduces the spread and shifting centroid of the metallicity distribution in each layer. As previously noted, the model predicts a slightly lower $[\alpha /\mathrm{Fe}]$ at the highest $[\mathrm{Fe}/{\rm{H}}]$, possibly a symptom of a ${t}_{\mathrm{sf}}$ that is too low, or of incorrect supernova yields, or incorrect APOGEE calibration at super-solar $[\mathrm{Fe}/{\rm{H}}]$. In the low-$| z|$ layer, the data show a hint of bimodality in $[\alpha /\mathrm{Fe}]$ at slightly sub-solar $[\mathrm{Fe}/{\rm{H}}]$, which is not predicted by the model. This hint of two distinct sequences could be a sign that radial migration has moved some low-α stars from the lower metallicity outer disk into the $R=3\mbox{--}5\,\mathrm{kpc}$ annulus (Schönrich & Binney 2009a), the complement to the explanation that H15 advance for the skew-positive shape of the midplane MDF at larger radii. Alternatively, the distinct low-α sequence could be connected to sharp changes in accretion history or star formation history (e.g., Chiappini et al. 1997; Few et al. 2014) not represented in the smooth models we consider here.

Zoom In Zoom Out Reset image size

Figure 8 displays the vertical distribution ${{dN}}_{* }/{dz}$ predicted by all nine models. For the central model, ${{dN}}_{* }/{dz}$ is well described (up to $| z| =4\,\mathrm{kpc}$) by the sum of two exponentials with scale-heights of 0.2 kpc and 0.6 kpc, with equal normalization at $| z| =1.036\,\mathrm{kpc}$. Although our model starts with ${z}_{h}=0.8\,\mathrm{kpc}$, few stars form at very early times, so the effective scale-height of the "thick disk" component is lower. We have arbitrarily normalized all of these distributions to ${{dN}}_{* }/{dz}=1\,{\mathrm{kpc}}^{-1}$ at z = 0. With the same midplane density, models that have longer t_c at a given ${t}_{\mathrm{sf}}$ produce more stars at high $| z|$, as expected. In the solar neighborhood, the Jurić et al. (2008) decomposition of the M-dwarf population implies thin- and thick-disk scale-heights of 0.3 kpc and 0.9 kpc, respectively, with equal normalizations at $| z| =0.95\,\mathrm{kpc}$. We do not know of any vertical density profile measurements at $R=3\mbox{--}5\,\mathrm{kpc}$ where our models apply, and the Jurić et al. (2008) sample is too local for reliable extrapolation to this Galactocentric radius. Such a measurement would be a valuable additional diagnostic for the scenario presented here. In particular, our choices of ${z}_{h}=0.8$ and 0.2 kpc for the initial and final scale-heights of our contracting disk, while loosely motivated by the Bird et al. (2013) cosmological simulation, are somewhat arbitrary. Additional constraints from a vertical profile measurement might allow these quantities to be treated as adjustable model parameters rather than fixed a priori.

Zoom In Zoom Out Reset image size

We have framed our model in terms of upside-down disk formation, with z_h(t) representing the scale-height of the star-forming gas layer over time. However, a model in which all stars formed in a thin layer and dynamical heating puffed up the older stellar populations would make identical predictions if it produced the same z_h(t). The principal argument that favors a contracting gas layer over stellar heating as the predominant mechanism governing disk thickness is that the z_h(t) of our most successful models is a plausible outcome of the upside-down scenario, with disk contraction occurring as the disk transitions from gas-dominated to having comparable gas and stellar mass. For ${\tau }_{* }=2.0\,\mathrm{Gyr}$, ${t}_{\mathrm{sf}}=3.0\,\mathrm{Gyr}$, and a recycling fraction r = 0.4, the ratio of stellar mass to gas mass is 0.5 at $t=2.5\,\mathrm{Gyr}$ and 1.2 at $t=4.5\,\mathrm{Gyr}$. In a heating scenario, these z_h(t) histories would require a heating mechanism that is sharply peaked at early times (see Figure 2), rather than the continuous, roughly power-law dependence on lookback time that is predicted by typical dynamical heating mechanisms (Hänninen & Flynn 2002). However, we have not explicitly demonstrated that a conventional heating model is inconsistent with the H15 data. Even in an upside-down scenario, dynamical heating is likely to influence the scale-height of populations at the present day (Bird et al. 2013, their Figure 19; J. Bird et al. 2017, in preparation).

We have not explicitly examined scenarios in which disk heating is a sudden event, as might result from a satellite merger. Given the drastic failure of our ${t}_{c}=1.0\,\mathrm{Gyr}$ models, however, we think it unlikely that a discrete heating scenario can reproduce the steady trend of the H15 abundance distributions with $| z|$.

It is possible that the average supernova yields we assume in our models are overestimated, if for example many massive stars collapse to form black holes or if the SN Ia rate is less than implied by Maoz et al. (2012). Because yields play an important role in determining equilibrium, we expect that lower assumed yields would affect which model is preferred by our analysis. A model with lower yields requires a lower value of η to achieve the same equilibrium abundance ratios as a model with higher yields (see Equations (8) and (9)). A lower η value implies a longer timescale for gas depletion and thus slower evolution toward equilibrium along a track in [O/Fe]–$[\mathrm{Fe}/{\rm{H}}]$. In order to quantify the impact of lower assumed yields, we have tried the simple experiment of halving ${m}_{{\rm{O}}}^{\mathrm{CC}}$, ${m}_{\mathrm{Fe}}^{\mathrm{CC}}$, and ${m}_{\mathrm{Fe}}^{\mathrm{Ia}}$, while adjusting η values for each ${t}_{\mathrm{sf}}$–t_c combination to give the same equilibrium ratios seen in our main set of models (i.e., slightly super-solar $[\mathrm{Fe}/{\rm{H}}]$). For ${t}_{\mathrm{sf}}=3\,\mathrm{Gyr}$, the low-yield model evolutionary track is nearly identical to the one generated by our central model. However, because lower yields result in slower evolution along the track, a collapse time of ${t}_{c}=4.5\,\mathrm{Gyr}$ would be preferred over ${t}_{c}=2.5\,\mathrm{Gyr}$ in the low-yield case.

Motivated by an "upside-down" scenario for thick disk formation, in which the height of the star-forming gas layer contracts as the stellar mass of the disk grows (Bournaud et al. 2009; Forbes et al. 2012; Bird et al. 2013), we have presented a simple model to explain the $| z|$-dependent $[\alpha /\mathrm{Fe}]\mbox{--}[\mathrm{Fe}/{\rm{H}}]$ distributions of inner disk stars ($3\,\mathrm{kpc}\,\lt R\lt 5\,\mathrm{kpc}$) measured by H15 from the APOGEE survey. The H15 inner disk stars appear to lie along a single evolutionary track, with most high-$| z|$ stars at sub-solar $[\mathrm{Fe}/{\rm{H}}]$ and $[\alpha /\mathrm{Fe}]=0.2\mbox{--}0.4$, most low-$| z|$ stars at super-solar $[\mathrm{Fe}/{\rm{H}}]$ and solar $[\alpha /\mathrm{Fe}]$, and mid-$| z|$ stars spread broadly over the $[\mathrm{Fe}/{\rm{H}}]$ range −0.5 to $+0.25$. Our model combines one-zone chemical evolution with a disk that contracts linearly from a scale-height of ${z}_{h}=0.8\,\mathrm{kpc}$ to ${z}_{h}=0.2\,\mathrm{kpc}$ over a timescale t_c, with ${z}_{h}=0.2\,\mathrm{kpc}$ thereafter. For specified assumptions about supernova yields and the SN Ia DTD, the shape of the model evolutionary track is governed mainly by the outflow mass-loading efficiency η, which sets the end-point in $[\mathrm{Fe}/{\rm{H}}]$ and the SFE timescale ${\tau }_{* }$, which sets the location of the knee where the track bends toward low $[\alpha /\mathrm{Fe}]$. The timescale ${t}_{\mathrm{sf}}$ of our linear-exponential star formation history, ${\dot{M}}_{* }(t)\propto {{te}}^{-t/{t}_{\mathrm{sf}}}$, has a small impact on the shape of the evolutionary track but a strong influence on the MDF. Together ${t}_{\mathrm{sf}}$ and the disk collapse timescale t_c determine the shape of the MDF in each $| z|$-layer. Our primary comparison of models to the H15 data appears in Figure 4.

The evolutionary track of the H15 sample is well matched by a model with ${\tau }_{* }=2.0\,\mathrm{Gyr}$, $\eta =1.5$, and ${t}_{\mathrm{sf}}=3.0\,\mathrm{Gyr}$. With ${t}_{c}=2.5\,\mathrm{Gyr}$, this model reproduces the observed MDF in each $| z|$-layer, and a model with ${t}_{c}=4.5\,\mathrm{Gyr}$ is almost as successful. With collapse timescale ${t}_{c}=1.0\,\mathrm{Gyr}$, however, the model predicts a strongly bimodal MDF at high $| z|$, contrary to observations. Models with sharply peaked early star formation, ${t}_{\mathrm{sf}}=1.0\,\mathrm{Gyr}$, predict a low-$| z|$ MDF that is much broader and more symmetric than the H15 measurement. Models with ${t}_{\mathrm{sf}}=6.0\,\mathrm{Gyr}$ yield good agreement with the low-$| z|$ MDF, but for any choice of t_c, they fail to match the observed MDF in the mid-$| z|$ and/or high-$| z|$ layer. Within our framework, therefore, the H15 data suffice to provide separate constraints on ${\tau }_{* }$, η, ${t}_{\mathrm{sf}}$, and t_c, at least at the factor-of-two level.

Our inferred z_h(t) could in principle be explained by either a contracting star-forming gas layer or age-dependent heating of stellar populations formed in a thin layer, or a combination thereof. A pure heating model seems unlikely to predict a z_h(t) that rises so sharply at large lookback times; for example, the scale-height of our ${t}_{c}=2.5\,\mathrm{Gyr}$ model triples between lookback times of $10\,\mathrm{Gyr}$ and $11.5\,\mathrm{Gyr}$ (or 8 and 11 Gyr for ${t}_{c}=4.5\,\mathrm{Gyr}$). Conversely, the inferred t_c corresponds to the interval over which the disk goes from being gas-dominated to having comparable gas and stellar mass fractions, a natural timescale for disk contraction in the upside-down scenario. The disk heating calculations of Aumer et al. (2016) predict a vertical velocity dispersion, and thus scale-height, that is roughly linear in age, which would correspond in our framework to a ${t}_{c}\approx 12.5\,\mathrm{Gyr}$ model. We have tried a variety of models with long t_c, and they all fail to match the H15 data.

The failure of our short-t_c models makes it unlikely that any scenario in which thick disk heating is a single discrete event can explain the continuous $| z|$-dependence of the MDF observed by H15. However, our model for z_h(t) is obviously idealized, and predictions of alternative models should be tested directly against the data. We do not claim that our model for the inner disk is unique, but it successfully describes the main characteristics of the H15 measurements with simple, physically plausible assumptions and a small number of adjustable parameters.

There are a number of directions for future progress in constraining and testing the scenario presented here. As illustrated in Figure 8, a direct measurement of ${{dN}}_{* }/{dz}$ for the inner disk ($R\approx 4\,\mathrm{kpc}$) would provide strong additional constraints, and perhaps allow the initial and final disk scale-heights to be treated as adjustable parameters without losing the power to constrain t_c and ${t}_{\mathrm{sf}}$. The DR13 APOGEE data set incorporates many improvements in abundance and stellar parameter determinations (SDSS Collaboration, 2016), and the APOGEE-2 survey of SDSS-IV (M. Blanton et al. 2017, in preparation) will eventually yield a substantially larger data set for studying the inner galaxy. More predictive semi-analytic (Forbes et al. 2012, 2014) or simulation-based (Bird et al. 2013; J. Bird et al. 2017, in preparation) models for the contracting gas layer could be combined with stellar heating prescriptions to produce more realistic z_h(t) templates and perhaps yield informative constraints on the physical parameters of such models. Full forward modeling, accounting for measurement errors and selection biases, and including nuisance parameters to represent uncertainties in supernova yields, would allow more rigorous statistical determination of parameter values and uncertainties. We include a simple analysis of the impact of APOGEE selection bias in Appendix B.

Stellar ages, inferred from asteroseismology (Pinsonneault et al. 2014; Stello et al. 2015) or chemical abundance signatures (Martig et al. 2016; Ness et al. 2016), would enable an especially direct test of our scenario. In the solar neighborhood, the age–metallicity relation exhibits a large dispersion (Edvardsson et al. 1993), a property naturally explained by models in which radial migration mixes stars formed at different Galactocentric radii (Wielen et al. 1996; Schönrich & Binney 2009a). We have focused on the inner disk because, in contrast to the solar neighborhood, its observed $[\alpha /\mathrm{Fe}]\mbox{--}[\mathrm{Fe}/{\rm{H}}]$ distribution resembles that predicted by a one-zone chemical evolution model with constant parameters. If our scenario is correct, then the age–metallicity relation of the inner disk should have low scatter and no dependence on $| z|$, and it should lie close to the relation predicted by our ${t}_{\mathrm{sf}}=3.0\,\mathrm{Gyr}$ models in Figure 3. While radial mixing by spiral arm perturbations is a symmetric process (Sellwood & Binney 2002), it may have much less impact on the inner disk because there are more stars available to migrate from small R to large R than vice versa (see H15, Figure 10). The horizontal ridge of stars visible in the low-$| z|$ $[\alpha /\mathrm{Fe}]\mbox{--}[\mathrm{Fe}/{\rm{H}}]$ distribution of H15 (bottom panel of Figure 1), while displaying only a hint of the bimodality in $[\alpha /\mathrm{Fe}]$ that becomes obvious at larger R, could be an indication that some of these stars have migrated inward from more distant birth sites. A complete test of our inner disk scenario should include a realistic allowance for radial migration from the outer disk.

As emphasized by Nidever et al. 2014 and H15, the "high-α locus" of stars in the $[\alpha /\mathrm{Fe}]\mbox{--}[\mathrm{Fe}/{\rm{H}}]$ plane is remarkably universal, with no detected dependence on R or $| z|$. This universality could indicate that our evolutionary scenario for the inner disk describes the Galaxy' early history over a much larger range of R, or that the high-α population throughout the disk is dominated by stars born at small R and scattered outward by radial migration. Assessing these possibilities requires embedding our model in a more complete framework that can account for the bimodality of the $[\alpha /\mathrm{Fe}]\mbox{--}[\mathrm{Fe}/{\rm{H}}]$ distribution observed over much of the Galactic disk. As surveys such as APOGEE, Gaia-ESO (Gilmore et al. 2012), LAMOST (Luo et al. 2015), and GALAH (De Silva et al. 2015) extend multi-element maps from the solar neighborhood to all regions of the Galaxy, the prospects for constructing, constraining, and testing comprehensive scenarios for the chemical evolution of the Milky Way grow dramatically brighter.

We thank Jennifer Johnson, Brett Andrews, Jonathan Bird, and Michael Hayden for valuable input at stages of this work. We thank Michael Hayden and Jon Holtzman for providing the distances used to define our observational comparison sample, and we acknowledge the immense efforts of the entire APOGEE team needed to produce this data set. This work was supported by NSF grant AST-60033987.

For our main analysis, we chose to adopt a linear-exponential form for the star formation (Equation (1)). An alternative and somewhat more conventional choice would be an exponential form, i.e.,

$$\begin{eqnarray}&&{\dot{M}}_{* }(t)\propto {e}^{-t/{t}_{\mathrm{sf}}}.\end{eqnarray} \tag{ 16 }$$

We consider a linear-exponential model to be more physically motivated, because it does not require starting with a large gas supply already in place at t = 0, and because it corresponds better to the star formation histories predicted in cosmological simulations (Simha et al. 2014). However, it is also interesting to see how our models perform if we adopt an exponential ${\dot{M}}_{* }(t)$.

Figure 9 repeats the comparison of Figure 4 for models with an exponential star formation history. Most of the trends seen in Figure 4 are repeated here. However, the low-$| z|$ MDF of the ${t}_{\mathrm{sf}}=3.0\,\mathrm{Gyr}$ exponential model fits the H15 data much less well than that of the corresponding linear-exponential model. The ${t}_{\mathrm{sf}}=6.0\,\mathrm{Gyr}$ model fares better, and with all layers considered, the parameter combination $({t}_{\mathrm{sf}},{t}_{c})=(6.0\,\mathrm{Gyr},2.5\,\mathrm{Gyr})$ is the most successful of the exponential models. It is somewhat less successful than our central linear-exponential model at reproducing the shapes of the mid-$| z|$ and high-$| z|$ MDFs. We conclude that our modeling provides moderate but not strong evidence for an inner disk star formation history that grows over time before turning downwards, rather than starting at a maximal rate.

Zoom In Zoom Out Reset image size

Our evolutionary models specify the number of stars born at each metallicity. In principle, the MDFs measured for evolved stars by APOGEE could differ significantly from these birth MDFs because the fraction of stars that are red giants of a given luminosity depends on the age of the stellar population. Because our one-zone models posit a one-to-one relation between age and metallicity, calculating the impact of APOGEE selection on the MDF measurement is fairly straightforward. Rather than simulate the field-by-field selection in detail, we apply weights to the observed luminosity distribution that reflect the fraction of populations of different ages in each magnitude bin. We make the following definitions.

• 1.  
  N_i is the number of observed stars in our sample (with the cuts defined in Section 2) in the ith bin of absolute H-band magnitude, $i=1\,...\,{n}_{L}$.
• 2.  
  M_j is the number of model stars formed in the jth bin of age, $j=1\,...\,{n}_{\mathrm{age}}$.
• 3.  
  F_ij is the fraction of all stars formed in the jth age bin that are in the ith bin of magnitude today, according to theoretical evolutionary tracks.
• 4.  
  A_kj is the fraction of stars in the jth age bin that are in the kth bin of $[\mathrm{Fe}/{\rm{H}}]$, $k=1\,...\,{n}_{\mathrm{FeH}}$. If our age bins corresponded to the timesteps of our model integration, then A_kj would equal one for a single value of k corresponding to the birth metallicity and zero for all other values of k. Because we choose broader age bins, A_kj is $\leqslant 1$ for a set of k and zero for most k.

With this notation, the normalized non-selection-adjusted MDFs, shown in Figure 4, may be written

$$\begin{eqnarray}&&{g}_{k}=\displaystyle \frac{{\displaystyle \sum }_{j}{A}_{{kj}}{M}_{j}}{{\displaystyle \sum }_{j}{M}_{j}}.\end{eqnarray} \tag{ 17 }$$

To compute the selection-adjusted MDFs, we first write the metallicity distribution of stars in the ith magnitude bin,

$$\begin{eqnarray}&&{f}_{k,i}=\displaystyle \frac{{\displaystyle \sum }_{j}{A}_{{kj}}{F}_{{ij}}{M}_{j}}{{\displaystyle \sum }_{j}{M}_{j}{F}_{{ij}}}.\end{eqnarray} \tag{ 18 }$$

Therefore, the selection-adjusted metallicity distribution across all magnitudes is

$$\begin{eqnarray}&&{\tilde{g}}_{k}=\displaystyle \frac{{\displaystyle \sum }_{i}{f}_{k,i}{N}_{i}}{{\displaystyle \sum }_{i}{N}_{i}}.\end{eqnarray} \tag{ 19 }$$

We compute ${f}_{k,i}$ for each of our models using luminosity functions from the Dartmouth Stellar Evolution Database (Dotter et al. 2008) to determine F_ij. M_j was obtained by binning the results of Equation (1), and similarly, binning stars formed in successive timesteps by age and metallicity gives A_kj. Figure 10 compares these selection-weighted MDFs with the purely theoretical ones considered in the body of this paper, showing that the corrections are small and do not affect our overall conclusions

Zoom In Zoom Out Reset image size