Inflow, Outflow, Yields, and Stellar Population Mixing in Chemical Evolution Models

The parameters of the fiducial simulation were chosen to broadly match the observed solar neighborhood trend in [O/Fe]–[Fe/H], and we use this set of parameters for all subsequent simulations unless stated otherwise. We adopt the net CCSN yields from Chieffi & Limongi (2004) and Limongi & Chieffi (2006), net asymptotic giant branch (AGB) star yields from Karakas (2010), and the W70 SN Ia yields from Iwamoto et al. (1999) (see Section 4 for more details). We normalized the abundances to the Lodders (2003) photospheric solar abundances.

The simulation starts with 2 × 10¹⁰ ${M}_{\odot }$ of primordial composition gas, and an additional 3.5 × 10¹¹ ${M}_{\odot }$ primordial composition gas flows into the galaxy with an exponentially declining time profile (τ = 6 Gyr). The final total mass of the simulation is 6.9 × 10¹⁰ ${M}_{\odot }$, with 3.2 × 10⁹ ${M}_{\odot }$ in gas, 5.2 × 10¹⁰ ${M}_{\odot }$ in stars, and 1.4 × 10¹⁰ ${M}_{\odot }$ in remnants (black holes, neutron stars, and WDs). Though we compare the model abundances to solar neighborhood data and match the model mass to the total mass of the Milky Way, the total mass of the model is essentially arbitrary (except for small stochastic supernova effects). All masses scale in proportion, so elemental abundances and gas fractions are independent of the total mass. Abundances across the Milky Way cannot be accurately modeled as a single zone due to different conditions at different radii. We present a suite of models in Section 5.2 as a first step toward capturing this diversity.

The SFR is set by a constant SFE of 10⁻⁹ yr⁻¹, converting ${M}_{\mathrm{gas}}(t)\cdot ({\rm{\Delta }}t/1$ Gyr) of the available gas mass into stars in a time interval ${\rm{\Delta }}t$. The cold interstellar medium (ISM) is ejected in an outflow at a rate of 2.5 × SFR. We adopt a Kroupa (2001) IMF from 0.1 to 100 ${M}_{\odot }$ and the stellar lifetime function of Renzini (1986, pp. 213–223) (see Equations (3) and (4) from Padovani & Matteucci 1993. The SN Ia DTD is an exponential with a timescale of 1.5 Gyr and a minimum delay time of 150 Myr. We list the parameters in Table 1 and explain the choice of parameters in more detail in Section 3.

Table 1.  Fiducial Model Parameters and Variations

Fiducial Parameters
Outflow mass loading factor	$\eta =2.5$
Star formation efficiency	SFE = 1.0 $\times {10}^{-9}$ yr−1
Exponential SN Ia DTD	${\tau }_{\mathrm{Ia}}=1.5$ Gyr, ${t}_{\min }=0.15$ Gyr
Exponential inflow history	${\tau }_{\inf }=6$ Gyr
IMF	Kroupa (2001) M = 0.1–100 ${M}_{\odot }$
Variations
Inflow timescale	Figure 3(a)
Inflow time history	Figures 3(b) and 4
Star formation efficiency	Figure 3(c)
Outflow mass loading	Figure 3(d)
Inflow metallicity	Figure 3(e)
IMF	Figure 3(f)
Minimum SN Ia delay time	Figure 5
Form of SN Ia DTD	Figures 6 and 7

Download table as:  ASCIITypeset image

To illustrate the importance of accretion and outflows, Figure 1 shows the evolution in [O/Fe]–[Fe/H], the MDF, and the [O/Fe] distribution function ([O/Fe]–DF) of the fiducial simulation, a closed box (no gas inflow or outflow) simulation, and an inflow-only simulation. For the purpose of comparison, we have adopted the fiducial simulation parameters for the closed box and inflow-only simulations, except for the gas flows themselves. Parameter adjustments could improve the agreement between these models and observations, though reproducing the observed track in full would require substantial changes to adopted supernova yields.

Zoom In Zoom Out Reset image size

The closed box simulation forms many of its stars early, as shown by the large high-α peak of its [O/Fe]–DF, and then steadily forms stars across a wide range of metallicities, producing a very broad MDF. It reaches an unreasonably high metallicity ([Fe/H] ≈ 1.5) because most of the gas has been turned into stars (${M}_{\star }^{\mathrm{final}}/{M}_{\mathrm{gas}}^{\mathrm{final}}=200$), leaving only a small amount of hydrogen to dilute the iron synthesized by CCSNe and SNe Ia.

The inflow-only simulation achieves a lower but still excessive metallicity of [Fe/H] ≈ 0.9 and finishes at a higher [O/Fe] (0.05 versus −0.05 for the closed box simulation), because inflow dilutes the metallicity by providing additional hydrogen and fuels more star formation and CCSN enrichment at late times relative to a closed box. This late time accretion results in the more prominent high metallicity and low-α peaks of the MDF and [O/Fe]–DF, respectively.

The fiducial simulation demonstrates the effect of outflows, which decrease the overall metallicity ([Fe/H] ≈ 0 at the end of the simulation by design) by removing metals and the gas that would have fueled more star formation and metal production. Outflows have the secondary effect of decreasing [O/Fe] because CCSN oxygen yields slightly increase with metallicity. The MDF and [O/Fe]–DF of the fiducial simulation are even more sharply peaked than the inflow-only simulation.

The fiducial simulation roughly reproduces the observed abundance trend in the data of Ramírez et al. (2013), though we caution that the detailed appearance of the observations depends on sample selection and on the methodology for abundance estimates. In the fiducial simulation and all of the other simulations that approximately reproduce observed abundance trends, the abundance evolution slows dramatically at late times, approaching an approximate equilibrium between sources and sinks of different elements. Below we refer to these approximately asymptotic states as equilibrium abundances. The emergence of equilibrium abundances and the timescales for reaching them can be understood analytically given some simplifying assumptions, as described by WAF.

Our goal in this paper is to understand how model predictions depend on model inputs, not to present a specific model of the Galaxy or the solar neighborhood. It is nonetheless worth explaining our choice of fiducial model parameters in the broader context of chemical evolution modeling, and especially our adoption of a high outflow mass loading factor $\eta =2.5$. Our choice of the stellar IMF is motivated by the observational evidence summarized by Kroupa (2001), and our choice of the SN Ia DTD is motivated by a combination of empirical and theoretical considerations as described in Section 3.6. With these choices, the predicted present-day ratio of CCSNe to SNe Ia in the fiducial model is 3.7, consistent with the (rather loose) empirical constraint of 4.3 ± 1.3 found by Li et al. (2011). Population averaged yields follow from the nucleosynthesis calculations cited in Section 2.1. Our adoption of a 6 Gyr e-folding timescale for gas infall is motivated by observational estimates of the star formation history (SFH) of the solar neighborhood, which suggest slowly declining star formation over the history of the Galaxy (Twarog 1980; Rocha-Pinto et al. 2000).

With these quantities set, the main adjustable parameters of the model are η and the SFE. As shown in Section 3 (see Figure 3), the value of the metallicity at late times is controlled mainly by η, and we choose η = 2.5 so that the fiducial model reaches [Fe/H] ≈ 0 at late times. For fixed η, the SFE controls the location of the knee in [O/Fe]–[Fe/H], and we select an SFE of ${10}^{-9}\,{\mathrm{yr}}^{-1}$ to reproduce the observed location at [Fe/H] ≈ −1.

The notion that strong outflows are required to achieve solar abundance in the ISM is common in the literature on the galaxy mass–metallicity relation (e.g., Finlator & Davé 2008; Peeples & Shankar 2011; Zahid et al. 2012), but it is not universal in the Galactic chemical evolution literature. Peeples et al. (2014) conclude that typical star-forming galaxies have ejected 75%–80% of the oxygen they have produced, and our fiducial model predicts 82% ejection. WAF show that in the instantaneous recycling approximation the oxygen abundance evolves to an equilibrium value ${Z}_{{\rm{O}},\mathrm{eq}}={m}_{{\rm{O}}}^{\mathrm{cc}}/(1+\eta -r)$ for a constant SFR (and slightly higher for a slowly declining SFR). Here ${m}_{{\rm{O}}}^{\mathrm{cc}}$ is the IMF-averaged oxygen yield from CCSNe and r is the mass recycling fraction, approximately 0.45 for a Kroupa IMF. Our IMF and yield assumptions imply ${m}_{{\rm{O}}}^{{cc}}\approx 0.017$ (i.e., $1.7\,{M}_{\odot }$ of oxygen produced per $100\,{M}_{\odot }$ of stars formed). The Lodders (2003) solar oxygen abundance, similar to that found in nearby B-stars (Przybilla et al. 2008), corresponds to ${Z}_{{\rm{O}}}=0.0056$, and reproducing it with ${m}_{{\rm{O}}}^{\mathrm{cc}}$ requires η = 2–3. We assume that all stars from 8 to 100 ${M}_{\odot }$ explode as CCSNe, and the oxygen yield could be significantly lower if the more massive stars (which produce most of the oxygen) instead collapse to black holes or if the wind loss from these stars is overpredicted. For example, if we cut out oxygen from stars with $M\gt 35\,{M}_{\odot }$ then ${m}_{{\rm{O}}}^{\mathrm{cc}}$ drops by a factor of two, reducing the required outflow efficiency to $\eta \approx 1$. However, one would also need to adjust the CCSN and SN Ia iron yields to reproduce solar iron abundance with this lower efficiency.

Many historical chemical evolution models have reproduced broad characteristics of the solar neighborhood without requiring outflows (e.g., Matteucci & Greggio 1986; Chiappini et al. 1997). Some of the earlier models effectively treated the population-averaged supernova yield as an adjustable parameter, concentrating on the form of the MDF rather than its normalization, and thus omitted the main constraint that drives us to substantial outflows. Another difference in many of these models is the use of a Salpeter (1955) IMF rather than a Kroupa (2001) or Chabrier (2003) IMF. The much larger population of low mass stars in a Salpeter IMF reduces the IMF-averaged supernova yield; for example, with the same supernova yields that we adopt, a Salpeter IMF gives ${m}_{{\rm{O}}}^{\mathrm{cc}}=0.011$ instead of 0.017, so the same oxygen abundance can be achieved with a lower η. Even with a Salpeter IMF, Matteucci & Greggio (1986) find that "the predicted solar iron abundance is higher by more than a factor of two than the observed one," which is precisely the mismatch that leads us to adopt strong outflows. The Chiappini et al. (1997) model incorporates both two episodes of infall and a decreased SFE during the second phase. It is unclear whether these changes can explain the lower abundances relative to a one-zone model with no outflow and fixed parameter; they may also be adopting lower yields than we do, but the descriptions are sufficiently different that a direct comparison is difficult. Radial gas flows (Spitoni & Matteucci 2011; Bilitewski & Schönrich 2012; Pezzulli & Fraternali 2016) are another possible mechanism for reducing metallicity by bringing in gas from the metal-poor outer disk. Nonetheless, in the context of one-zone models with an IMF and supernova yields similar to those adopted here, it is clear that strong outflows are a necessary ingredient for reproducing solar abundances.

For simplicity, we have restricted our analysis to models in which the outflow efficiency and SFE are constant in time, though either could plausibly evolve as the disk potential grows and the gas fraction declines. The high SFE of our fiducial model, chosen to reproduce the location of the [α/Fe]–[Fe/H] knee, leads to a low gas fraction at late times, ${M}_{\mathrm{gas}}/{M}_{* }\approx 0.05$. This is considerably lower than the recent estimate ${M}_{\mathrm{gas}}/{M}_{* }\approx 0.4$ for the solar neighborhood by McKee et al. (2015). For a specified SFH, the gas fraction is inversely proportional to the SFE, so this mismatch likely indicates a substantially decreased SFE at late times, relative to the high efficiency at early times when [α/Fe] is just beginning its decline toward solar values.

Two other solar-neighborhood constraints we do not attempt to match are the age–metallicity relation (Twarog 1980; Edvardsson et al. 1993; Casagrande et al. 2011) and the MDF (Pagel & Patchett 1975; Nordström et al. 2004; Casagrande et al. 2011). Both of these quantities are likely to be strongly shaped by radial migration of stars, as argued by Schönrich & Binney (2009a) and Hayden et al. (2015) building on suggestions by Wielen (1977) and Sellwood & Binney (2002). One-zone models with constant parameters predict a one-to-one age–metallicity relation, while the observed relation shows large scatter (Edvardsson et al. 1993). As shown below, these models generically predict negatively skewed, sharply peaked MDFs like those seen in the inner Galaxy (Hayden et al. 2015), not the broad, roughly Gaussian (in [Fe/H]) form seen in the solar neighborhood. One-zone models provide a useful starting point for inferring the impact of radial migration, which mixes populations formed with different enrichment histories, but a single one-zone model cannot explain these aspects of solar neighborhood chemical evolution.

Since we will be using [O/Fe]–[Fe/H] as the main diagnostic plot, it is helpful to understand the enrichment sources responsible for producing oxygen and iron. Figure 2 shows the mass of gas-phase oxygen and iron as a function of metallicity (stellar abundances are taken to be the gas-phase abundances when a stellar generation is born). For ease of interpretation we show the results from a constant SFR simulation because its iron mass is linear with [Fe/H]; this simulation follows a very similar locus in [O/Fe]–[Fe/H] to the fiducial simulation, though its final [Fe/H] is lower and [O/Fe] is higher (see Figure 3(b)). The solid line shows the total gas-phase mass of each element. The dotted, short-dashed, and long-dashed lines indicate the contributions from CCSNe, AGB stars, and SNe Ia, respectively. The ordinate axes have been offset such that the total oxygen and iron lines cross at solar [O/Fe].

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

One can metaphorically regard CCSN, SN Ia, and AGB stars as tributaries that enrich the "river" of the main gas reservoir, which is itself depleted by star formation and outflow. Far "upstream," at early times and low [Fe/H], CCSNe are the only significant channel and the ratio [O/Fe] ≈ 0.4 reflects their high-α yields. As SNe Ia become important, contributing iron but minimal oxygen, the [O/Fe] ratio declines toward solar values. In the constant SFR simulation, the transition from CCSN to SN Ia dominating iron production occurs at t = 2.7 Gyr and [Fe/H] ≈ −0.17. CCSNe remain the dominant source of oxygen production at all times, with a much smaller contribution from AGB stars. Since these lines include recycled material, AGB stars "contribute" iron even though they do not directly synthesize it because they return iron incorporated at birth. For similar plots of the full suite of elements and more in-depth discussion of multi-element abundances see Figure 10 and Section 4.2.

Closed box chemical evolution models generically suffer from the G-dwarf problem (van den Bergh 1962), namely that they produce MDFs for G-dwarfs (whose lifetimes are comparable to the age of the Galactic disk) with too many metal-poor stars relative to metal-rich stars. Inflows allow galaxies to form more stars later in their lifetimes when the cold ISM is metal-rich, increasing the fraction of metal-rich stars. Continuing infall is, of course, the natural expectation in analytic or numerical models of cosmological galaxy formation, with star formation typically tracking gas accretion after a moderate delay (e.g., Katz et al. 1996).

To represent a range of inflow histories, we vary the inflow timescale ${\tau }_{1}$ assuming an exponential inflow rate:

$$\begin{eqnarray}&&{\dot{M}}_{\mathrm{in}}=\displaystyle \frac{{M}_{1}}{{\tau }_{1}}{e}^{-t/{\tau }_{1}},\end{eqnarray} \tag{ 1 }$$

where M₁ sets the normalization of the inflow rate. Figure 3(a) shows the effect of varying the inflow timescale by a factor of two relative to the ${\tau }_{1}=6$ Gyr of the fiducial simulation. All three simulations follow very similar tracks, but the equilibrium abundances, [O/Fe]–DFs, and MDFs are distinct. The ${\tau }_{1}=3$ Gyr simulation has a lower equilibrium [O/Fe] and a higher equilibrium [Fe/H]. It experiences less late-time accretion, star formation, and CCSN enrichment to counteract SN Ia Fe production, which continues to drive down [O/Fe] and increase [Fe/H]. Its [O/Fe]–DF and MDF are broader because it does not asymptote to an equilibrium abundance as quickly as the fiducial simulation. Conversely, the ${\tau }_{1}=12$ Gyr simulation reaches a higher [O/Fe] and lower [Fe/H] at the end of the simulation and has more peaked [O/Fe]–DF and MDF. Overall, the inflow timescale dictates how far along the track the simulation goes until it achieves equilibrium abundance ratios at late times. The faster a simulation reaches the equilibrium abundance, the more peaked its [O/Fe]–DF and MDF will be.

We experimented with several functional forms of the inflow history with different ratios of early to late time accretion, including a double exponential,

$$\begin{eqnarray}&&{\dot{M}}_{\mathrm{in}}=\displaystyle \frac{{M}_{1}}{{\tau }_{1}}{e}^{-t/{\tau }_{1}}+\displaystyle \frac{{M}_{2}}{{\tau }_{2}}{e}^{-t/{\tau }_{2}},\end{eqnarray} \tag{ 2 }$$

and a linear–exponential product,

$$\begin{eqnarray}&&{\dot{M}}_{\mathrm{in}}=\left(\displaystyle \frac{{M}_{1}}{{\tau }_{1}}\right)\left(\displaystyle \frac{t}{{\tau }_{1}}\right){e}^{-t/{\tau }_{1}}.\end{eqnarray} \tag{ 3 }$$

We also ran a constant SFR simulation, where the inflow rate was set to maintain a constant gas mass at all times. This simplified SFH helps illustrate the differences between simulations with different yields (see Figure 9) and the relative contribution of CCSN, SN Ia, and AGB stars to enrichment (see Figures 2 and 10).

Figure 4 shows the SFR as a function of time for the fiducial, double exponential, linear–exponential, and the constant SFR simulations. The double exponential simulation represents a scenario with early rapid accretion (M₁ = 10¹⁰ ${M}_{\odot }$, ${\tau }_{1}=300$ Myr) followed by steady accretion on long timescales (M₂ = 6 × 10¹¹ ${M}_{\odot }$, ${\tau }_{2}=14$ Gyr), where the timescales were adopted from Schönrich & Binney (2009a). The linear–exponential simulation is motivated by cosmological hydrodynamical simulations that find inflow rates that increase at early times, reach a maximum at a characteristic timescale, and then decline at late times (Simha et al. 2014), so it starts with no initial gas. The timescale of the linear–exponential simulation (${\tau }_{1}=3.5$ Gyr) was chosen so that the simulation reached solar metallicity and oxygen abundance. The normalizations of the double exponential, linear–exponential, and the constant SFR simulations were chosen to approximately match the same total amount of inflow and final stellar mass as the fiducial simulation.

Zoom In Zoom Out Reset image size

Figure 3(b) shows the mean tracks in [O/Fe]–[Fe/H], the [O/Fe]–DFs, and MDFs for these inflow histories. The double exponential simulation enriches slightly faster than the fiducial simulation for the first 2.5 Gyr, but its higher late-time inflow rate results in a higher equilibrium [O/Fe], a lower equilibrium [Fe/H], and more peaked [O/Fe]–DF and MDF. The linear–exponential simulation starts at a lower metallicity ([Fe/H] = −2.0 at t = 30 Myr) than the fiducial simulation because it does not start with any gas. Its lower early-time SFR shifts the knee to lower metallicity and produces fewer stars on the high [O/Fe] plateau. It also broadens the [O/Fe]–DF and MDF, especially toward the higher [O/Fe] and lower [Fe/H] side of the peak. However, its equilibrium [O/Fe] and [Fe/H] match those of the fiducial simulation. Since the constant SFR simulation effectively has an infinite inflow timescale, it behaves similarly to the τ = 12 Gyr inflow timescale simulation. The constant SFR simulation follows nearly the same track as the fiducial simulation, but its equilibrium [O/Fe] is higher, equilibrium [Fe/H] is lower, and its [O/Fe]–DF and MDF are more peaked due to the additional dilution and CCSN enrichment at late times from its higher late inflow rate.

The total mass of the simulation does not affect the mean trend in [O/Fe]–[Fe/H], as long as the initial mass and the mass inflow rate are scaled together. However, larger masses decrease the effect of stochastic fluctuations from incomplete sampling of the IMF or the SN Ia DTD. For a fixed inflow rate, the initial gas mass of the model regulates the speed of enrichment but not the final metallicity. In this case, a smaller initial gas mass extends the plateau in [O/Fe]–[Fe/H] to lower metallicities. A large initial gas mass produces a galaxy that enriches quickly then asymptotes to a constant metallicity, whereas a small initial gas mass results in a smoother, more continuous increase in metallicity.

A common prescription for SFR is the Kennicutt–Schmidt law (${\dot{{\rm{\Sigma }}}}_{\star }\propto {{\rm{\Sigma }}}_{\mathrm{gas}}^{N};$ Schmidt 1959; Kennicutt 1998). Kennicutt (1998) found a best-fit power-law index of N = 1.4 across a wide range in global galaxy SFRs, which implies a higher SFE in denser regions. More recent studies by Leroy et al. (2008) and Bigiel et al. (2008) determined a value of N = 1.0 for molecular-dominated gas in spatially resolved observations of local galaxies, which implies a constant SFE for gas that is sufficiently dense. Scoville et al. (2016) also found a nearly linear slope (N = 0.9) for high redshift (z = 1–6) star-forming galaxies using the dust continuum to measure the total gas mass. Likewise, our model assumes a constant SFE (i.e., N = 1.0):

$$\begin{eqnarray}&&\mathrm{SFR}=\displaystyle \frac{1}{{t}_{\mathrm{gas}}}\,{M}_{\mathrm{gas}}=\mathrm{SFE}\cdot {M}_{\mathrm{gas}},\end{eqnarray} \tag{ 4 }$$

where t_gas is the gas depletion timescale and $\mathrm{SFE}\equiv {t}_{\mathrm{gas}}^{-1}$ is the SFE. This parameterization is the natural choice for models with a constant SFR or if the gas is predominantly molecular, but for models with a declining gas supply and a non-negligible fraction of atomic gas, a declining SFE is arguably a more observationally motivated assumption. That being said, we tested simulations with a Kennicutt–Schmidt star formation law (N = 1.4), and their evolution was not significantly different from the fiducial simulation. We note that Portinari & Chiosi (1999) found that a Kennicutt–Schmidt law with either N = 1.0 or 1.5 cannot simultaneously produce the radial metallicity gradient and gas profile, though radial gas flows can alleviate this issue (Portinari & Chiosi 2000). As for the normalization of the SFE, Leroy et al. (2008) and Bigiel et al. (2008) measured the SFE for molecular gas to be 5 × 10⁻¹⁰ yr⁻¹, which is a factor of two lower than the SFE required in the fiducial simulation to match the knee in [O/Fe]–[Fe/H]. However, this comparison is dependent on the adopted minimum SN Ia delay time (see Figure 5), so it is inadvisable to draw strong conclusions about the agreement between the overall normalization of the SFE in observations and the model.

Zoom In Zoom Out Reset image size

Figure 3(c) shows factor-of-three variations around the SFE of the fiducial simulation (10⁻⁹ yr⁻¹). Increasing the SFE results in a knee at a higher [Fe/H] because more CCSNe produce more iron and more star formation consumes more hydrogen (see also Matteucci & Brocato 1990). Further increases in the SFE result in diminishing increases in the [Fe/H] of the knee because the simulation reaches a temporary equilibrium abundance. Increasing the SFE also increases the relative number of stars on the high [O/Fe] plateau. The final equilibrium [O/Fe] and [Fe/H] are mostly insensitive to SFE, though SFE sets the speed of convergence, hence the width of the peaks of the [O/Fe]–DF and MDF become narrower with increasing SFE. To simplify the SFH, we also performed constant SFR simulations with the same three SFEs and find similar general characteristics for the SFE variations (modulo the effects of using a constant SFR mentioned above).

The energy and momentum injection from massive stars can launch large-scale galactic winds (e.g., Veilleux et al. 2005). Measurements of the mass outflow rates of these winds, while challenging, point to values that are a few times the SFR (e.g., Martin 1999; Heckman et al. 2000). We parameterize the outflow rate (${\dot{M}}_{\mathrm{outflow}}$) to be proportional to the SFR,

$$\begin{eqnarray}&&{\dot{M}}_{\mathrm{outflow}}={\eta }_{\mathrm{wind}}\mathrm{SFR},\end{eqnarray} \tag{ 5 }$$

where ${\eta }_{\mathrm{wind}}$ is the mass-loading factor of the wind. We assume that the outflowing material has the same abundance pattern as the cold ISM, which would be expected if the mass in the outflow is dominated by entrained cold gas as opposed to the hot gas ejected by SNe. We performed simulations that ejected a fraction of the stellar yields without any cold ISM leaving the Galaxy, which required 75% of the stellar yields to be ejected for the simulation to have a equilibrium [O/Fe] and [Fe/H] at solar values. However, the [O/Fe] knee of this simulation occurred at [Fe/H] = −1.5, which is significantly lower than the observed metallicity of the knee. Metallicity-enhanced winds can also spoil the otherwise good agreement of model predictions with the observed ISM oxygen abundance and deuterium-to-hydrogen ratio (Weinberg 2016).

Figure 3(d) shows the mean track in [O/Fe]–[Fe/H], the [O/Fe]–DF, and the MDF for factor of two variations in the mass-loading factor around the fiducial value of ${\eta }_{\mathrm{wind}}=2.5$, which produces a track that ends at solar metallicity and oxygen abundance. Changing ${\eta }_{\mathrm{wind}}$ strongly affects the shape of the track after the knee at [Fe/H] ≈ −0.9, leading to very different equilibrium abundances, especially [Fe/H]. Increasing the mass-loading factor removes more metal-rich gas, which is replaced by inflowing pristine gas, thus decreasing the equilibrium metallicity. The lower equilibrium metallicity results in a lower equilibrium [O/Fe] because the CCSN oxygen yields slightly increase with metallicity. The [O/Fe]–DF and MDF of the three simulations are offset because of their different equilibrium abundances, but the shape of the peaks remain similar. A higher mass-loading factor also induces greater variance in the abundances because it decreases the size of the gas reservoir and boosts the impact of stochastic fluctuations, but this effect only becomes important at much lower total masses than shown in Figure 3(d).

Gas that accretes onto galaxies is a mix of pristine or low metallicity intergalactic medium (IGM) gas (Z ≈ 10⁻³ ${Z}_{\odot }$; Songaila & Cowie 1996) and enriched material previously ejected from galaxies (e.g., Marasco et al. 2012; Ford et al. 2014). The composition of accreting gas is difficult to determine observationally, so the relative contribution of new versus recycled gas is unknown. We performed simulations where the overall metallicity of the inflowing gas was

• 1.  
  Z = 0 (fiducial value used to simulate inflow of pristine gas),
• 2.  
  Z = 10⁻² ${Z}_{\odot }$ (to simulate inflow from an enriched IGM), and
• 3.  
  Z = 10⁻¹ ${Z}_{\odot }$ (to simulate a mix of pristine gas and recycled gas),

where ${Z}_{\odot }$ = 0.02 (this value is higher than the sum of the solar abundances of individual elemental abundances of Lodders 2003 that we have adopted, though this difference does not affect our qualitative findings). The composition of the inflow was the relative abundances of elements in the ISM at the previous time step, which was then scaled to an overall metallicity. We also experimented with α-enhanced and scaled-solar abundance patterns but chose not to use them because the α-enhanced abundance pattern prevented the gas from reaching the solar [α/Fe] value and the scaled-solar abundance pattern created very low metallicity stars with solar [α/Fe] abundances.

Figure 3(e) shows the simulations with enriched inflow. The Z = 10⁻² ${Z}_{\odot }$ simulation is nearly identical to the fiducial (Z = 0) simulation, indicating that enriched inflow from the IGM has little effect on chemical evolution (though it may be important at extremely low metallicities; see Brook et al. 2014). The Z = 10⁻¹ ${Z}_{\odot }$ simulation has a higher [Fe/H] and a higher [O/Fe] (above the knee) at fixed time than the fiducial simulation, so the knee occurs at a higher [Fe/H]. Its [O/Fe]–DF and MDF are peaked at higher values, though the shape of these distribution functions is similar to the fiducial simulation. We conclude that enriched inflow only has a major effect on the track in [O/Fe] versus [Fe/H] if the accreted gas has a metallicity of $Z\gt {10}^{-1}$ ${Z}_{\odot }$.

In this experiment, we have only considered constant inflow metallicities that stay fixed throughout the simulations. In reality, the inflow metallicity likely increases with time and would be lower than the ISM metallicity at early times. This issue is significantly mitigated by the dominance of the large reservoir of pristine gas relative to the early infalling gas, so the infall metallicity has a small impact on the ISM metallicity provided that the infall metallicity is low compared to the eventual equilibrium metallicity. This finding is confirmed in the analytic model of WAF (see their Section 5.2). Cosmological simulations (e.g., Kereš et al. 2009) find that the transition from pristine cold flows to enriched hot mode accretion occurs around z ∼ 1, when the abundances are already approaching equilibrium values. Thus, we expect that a time-dependent inflow metallicity will have a negligible effect in moderate and high metallicity environments, though it may be important in the low metallicity outer regions of galactic disks.

The slope and mass range of the stellar IMF sets the relative number of high and low mass stars that form in a stellar population. The main effect of the IMF is that it governs the relative mass locked up in low mass stars and remnants versus the element production by high and intermediate mass stars. Additionally, the high mass slope and cutoff can have a secondary effect on elements whose CCSN yields are highly mass dependent, like oxygen. Côté et al. (2016a) found that the high mass slope of the IMF is a major source of uncertainty in chemical evolution models.

Historically, the Salpeter (1955) IMF, with a slope of α = 2.35 and a mass range from 0.1 to 100 ${M}_{\odot }$, has provided a useful reference point for IMFs. However, the more recent Kroupa (2001) IMF, whose slope is α = (1.3, 2.3) from (0.1–0.5, 0.5–100 ${M}_{\odot }$) does a better job of describing the solar neighborhood, so we adopt it as the fiducial IMF. We have assumed that the IMF is constant in time, as suggested by the chemical evolution study of Romano et al. (2005).

In Figure 3(f), we show the mean track in [O/Fe]–[Fe/H], the [O/Fe]–DF, and the MDF for four IMFs:

• 1.   
  Kroupa IMF from 0.1 to 100 ${M}_{\odot }$ (fiducial),
• 2.   
  flat (α = 2.1) IMF from 0.1 to 100 ${M}_{\odot }$,
• 3.   
  Kroupa IMF from 0.1 to 50 ${M}_{\odot }$, and
• 4.  
  Salpeter IMF from 0.1 to 100 ${M}_{\odot }$.

Relative to the fiducial IMF, the flat IMF simulation always has a higher [O/Fe]. It also has a slightly higher [Fe/H], especially at early times ($t\lt 2$ Gyr). The Kroupa 0.1–50 ${M}_{\odot }$ IMF simulation has a lower [O/Fe] but a similar [Fe/H] compared to the fiducial simulation. These changes are due to the oxygen yield of CCSN increasing strongly with mass but iron yield staying constant.

The Salpeter 0.1–100 ${M}_{\odot }$ IMF simulation has a nearly identical [O/Fe] value of the plateau because of the similar high mass slopes of the Kroupa and Salpeter IMFs and the lack of SN Ia iron production at these early time steps. However, the larger number of low mass stars that are produced in the Salpeter IMF simulation results in a lower [Fe/H] because more gas goes into long-lived low mass stars as opposed to CCSNe. The [O/Fe]–DFs and MDFs for the IMF variation simulations have similar shapes, though the locations of their peaks are offset for the reasons given above.

Although Figure 3(f) shows that chemical evolution tracks are sensitive diagnostics of the IMF, these differences are partly degenerate with SN yields (see Mollá et al. 2015). The mass cut for CCSNe, which sets the division between mass ejected in the explosion and mass that falls back on the neutron star, has a large impact. In particular, the iron yield is governed by the mass cut because the material near the mass cut consists mostly of iron-peak elements, whereas oxygen is produced further out in the star so its yield is insensitive to the mass cut. Increasing the mass cut to produce 0.05 ${M}_{\odot }$ of ejected ⁵⁶Ni instead of 0.1 ${M}_{\odot }$, increases the [O/Fe] level of the plateau, decreases the [Fe/H] of the knee, and generally decreases [Fe/H] at all times. It also produces slightly broader [O/Fe]–DFs and MDFs. We discuss mass cut effects further in Section 4.1.1.

We note that the [O/Fe]–[Fe/H] tracks for the various IMFs assume the outflow mass loading parameter for the Kroupa IMF. The Kroupa IMF is more bottom light/top heavy than the alternative IMFs that we consider here and the IMFs others have considered elsewhere (e.g., Scalo 1986; Kroupa & Weidner 2003). Recalibrating the model for one of these IMFs would lead to a smaller outflow mass loading parameter, but at least in the case of the Salpeter IMF from 0.1 to 100 ${M}_{\odot }$, the model still requires outflows to produce the knee at [Fe/H] = −1 and reach solar abundances.

The SN Ia history cannot be predicted robustly due to uncertainty in the delay time between the birth of a stellar population and the SN Ia that it will produce. Our understanding of the SN Ia DTD is complicated by a few factors: (1) the evolutionary state of the progenitor systems—double degenerate (WD + WD) or single degenerate (WD + red giant, WD + main sequence star, or WD + helium star)—remains unknown; (2) binary evolution, especially the common envelope phase, is difficult to model; (3) dynamical effects in triple star systems could dramatically hasten the merger of two WDs relative to naive expectations from binary evolution (Thompson 2011; Katz & Dong 2012).

Because of these uncertainties, we consider variations in the minimum delay time and multiple functional forms of the DTD. The minimum delay time sets the timescale for SNe Ia to become significant contributors of iron, which produces the knee in [O/Fe]–[Fe/H]. The different functional forms of the DTD govern the SN Ia rate on longer timescales.

In the fiducial simulation, we adopt an exponential form of the SN Ia DTD as in Schönrich & Binney (2009a), where the mass in the WD reservoir produced by a stellar population is depleted at a constant fractional rate after a minimum delay time. The SN Ia DTD (in units of SN Ia yr⁻¹ ${M}_{\odot }$⁻¹) is

$$\begin{eqnarray}&&{f}_{\mathrm{SN}\mathrm{Ia}}(t)=K\,{\tau }_{\mathrm{SN}\,\mathrm{Ia}}^{-1}\,{e}^{-(t-{t}_{\mathrm{SN}\mathrm{Ia}}^{\min })/{\tau }_{\mathrm{SN}\mathrm{Ia}}},\end{eqnarray} \tag{ 6 }$$

where ${t}_{\mathrm{SN}\,\mathrm{Ia}}^{\min }$ is the minimum delay time for a SN Ia and ${\tau }_{\mathrm{SN}\mathrm{Ia}}$ is the timescale for SNe Ia, which are set to 150 Myr and 1.5 Gyr from the birth of a stellar population, respectively, in the fiducial simulation (also adopted from Schönrich & Binney 2009a, who calibrated these values to reproduce oygen abundances). The constant K was chosen to match the observed time-integrated (from 40 Myr to 10 Gyr) number of SNe Ia per stellar mass formed ${N}_{\mathrm{Ia}}/{M}_{\star }=2.2\times {10}^{-3}\,{M}_{\odot }^{-1}$ from Maoz & Mannucci (2012). These authors adopted the Bell et al. (2003) IMF for their stellar mass estimates, which is similar enough to the Kroupa IMF for our purposes, but one would have to recalibrate the value of K if using a significantly different IMF.

Figure 5 shows the mean tracks in [O/Fe]–[Fe/H], the [O/Fe]–DFs, and the MDFs for an exponential SN Ia DTD with three different minimum delay times: 40, 150 (fiducial value), and 300 Myr. The fiducial simulation turns over at [Fe/H] ∼ −0.9. Increasing the minimum delay time increases the [Fe/H] of the knee and makes the turnover sharper. After about 2 Gyr, the [O/Fe]–[Fe/H] tracks of the three simulations are very similar. The MDFs of the three simulations are almost identical, though the relative number of stars on the high [O/Fe] plateau increases with the minimum delay time, as can be seen in the high [O/Fe] peak of the [O/Fe]–DF. The data do not rule out any of these possible delay times.

We tested three additional functional forms of the SN Ia DTD: a power-law DTD, a prompt + delayed DTD, and a theoretical population synthesis DTD, which are compared in Figure 6. The power-law DTD is motivated by recent measurements of the SN Ia rates and the SFHs of galaxies that have placed interesting observational constraints on the SN Ia DTD using a variety of strategies (Totani et al. 2008; Maoz & Badenes 2010; Maoz et al. 2010, 2011, 2012; Graur et al. 2011; Maoz & Mannucci 2012). Encouragingly, the results from different techniques agree that a ${t}^{-1}$ power law DTD provides a good fit to the data,

$$\begin{eqnarray}&&{f}_{\mathrm{SN}\mathrm{Ia}}(t)\propto {\left(\displaystyle \frac{t}{\mathrm{yr}}\right)}^{-1}.\end{eqnarray} \tag{ 7 }$$

The evidence for a ${t}^{-1}$ power-law DTD is strongest for t = 1–10 Gyr (Maoz & Mannucci 2012), but the data are consistent with a ${t}^{-1}$ power law continuing to shorter delay times with a minimum delay time of 40 Myr from the birth of a stellar population, which is approximately the lifetime of the most massive star that will produce a WD. We also consider an empirical prompt + delayed DTD based on observations that the SN Ia rate is correlated with SFR tracers (Scannapieco & Bildsten 2005; Mannucci et al. 2006). In this scenario, the "prompt" SNe Ia are associated with young stellar populations, so their rate scales with the SFR. The "delayed" SNe Ia come from old stellar populations, so their rate depends on stellar mass. The SN Ia rate is

$$\begin{eqnarray}&&\displaystyle \frac{{R}_{\mathrm{SN}\mathrm{Ia}}(t)}{{\mathrm{yr}}^{-1}}=A\left(\displaystyle \frac{{M}_{\star }}{{M}_{\odot }}\right)+B\left(\displaystyle \frac{\mathrm{SFR}}{{M}_{\odot }\,{\mathrm{yr}}^{-1}}\right),\end{eqnarray} \tag{ 8 }$$

where A = 4.4 × 10⁻¹⁴ and B = 2.6 × 10⁻³ are the values from Scannapieco & Bildsten (2005) converted to the relevant units, and we have implemented a 40 Myr minimum delay time for WD formation. Maoz & Mannucci (2012) point out that this scenario is similar to coarse time sampling of a ${t}^{-1}$ power-law DTD, and the demarcation between the prompt and delayed components is somewhat arbitrary. Finally, we implement the theoretical DTD from Greggio (2005) (described in their Section 3.1) based on stellar population synthesis models of single degenerate systems that reach the Chandrasekhar limit.

Zoom In Zoom Out Reset image size

Figure 6(a) shows the DTDs for the exponential, power law, and population synthesis DTDs for a stellar population born at t = 0, normalized to the observed time-integrated (from 40 Myr to 10 Gyr) SN Ia rate of ${N}_{\mathrm{Ia}}/{M}_{\star }=2.2\times {10}^{-3}\,{M}_{\odot }^{-1}$ (Maoz & Mannucci 2012; though see Côté et al. 2016a for a discussion of the significant element-dependent impact caused by uncertainty in this normalization). The power-law DTD produces more SN Ia at early times (≲200 Myr) and late times (≳4 Gyr) than the exponential DTD but fewer in between. The population synthesis DTD has a roughly similar shape to the exponential DTD, except shifted to shorter delay times. The breaks in the population synthesis DTD are due to the requirement of a carbon–oxygen WD and that the system must achieve the Chandrasekhar mass.

Figure 6(b) shows the number of SNe Ia per 30 Myr time step for simulations (with continuous star formation) with different DTDs but the same SFH. The numbers of SNe Ia per time step in the exponential, power law, and population synthesis DTD simulations increase up to a peak at 1.5–3 Gyr followed by a decline at late times. The prompt + delayed scenario produces a very high initial rate because it is tied to the SFR, which is high at early times. We adopted the normalization of Scannapieco & Bildsten (2005), though these numbers were determined for galaxies in equilibrium. The large number of very early SNe Ia has a dramatic effect in the [O/Fe]–[Fe/H] diagram. After 12 Gyr, the prompt + delayed scenario produces 30% more SNe Ia than the fiducial DTD due to the normalization. The exponential DTD with ${t}_{\min }=40$ Myr and the population synthesis scenario have nearly the same number of SNe Ia as the fiducial DTD after 12 Gyr. The power law DTD has 5% fewer SNe Ia after 12 Gyr than the fiducial DTD even though they were normalized to the same value because of the different relative number of short versus long delay SNe Ia. Another way to understand this point is that there are more SeN Ia with delays of 4–10 Gyr for the power law DTD than the fiducial DTD (see Figure 6(a)). For populations between 4 and 10 Gyr, more of their SNe Ia will have exploded by the end of the simulation for the fiducial DTD than the power-law DTD.

Figure 7 shows the mean tracks in [O/Fe]–[Fe/H], the [O/Fe]–DFs, and the MDFs of simulations with various SN Ia DTDs. We also include an exponential DTD with ${t}_{\min }^{\mathrm{SN}\,\mathrm{Ia}}$ = 40 Myr (blue dashed line) to more closely match the other DTDs. The power-law DTD simulation turns over at a lower metallicity than the exponential DTD, but it reaches the same [Fe/H] and [O/Fe] at about 2 Gyr. It then continues along the same track at a slower pace but ends at a higher equilibrium [Fe/H] and lower [O/Fe] than the exponential DTD. Its [O/Fe]–DF and MDF are significantly broader than those of the fiducial simulation, especially the [O/Fe]–DF, and do not show a sharp cutoff. The power-law DTD does not reach an equilibrium metallicity or [O/Fe] because of the large number of SNe Ia at late times (see also Kubryk et al. 2015a). The population synthesis DTD has a similar track to the power law DTD over the first 500 Myr. It reaches the same track as the exponential DTD but finishes with a lower equilibrium [Fe/H] and a higher [O/Fe] (near the 4 Gyr symbol). Its [O/Fe]–DF and MDF are more peaked than the fiducial simulation's. The exponential, power-law, and population synthesis scenario DTD are all consistent with the [O/Fe]–[Fe/H] data, especially if other parameters, such as the inflow timescale or the SN Ia normalization, are allowed to vary to compensate for the different equilibrium [Fe/H] and [O/Fe]. On the other hand, the prompt + delayed DTD creates a sudden drop in [O/Fe] before leveling off, a feature that is not present in the data. This behavior is caused by the large number of prompt SNe Ia (see Figure 6(b)). Its [O/Fe]–DF and MDF are strongly peaked, and its [O/Fe]–DF has a sharp left edge, which is unique among the models and inconsistent with the data. Matteucci et al. (2009) found that using an alternative parameterization of the prompt + delayed DTD from Mannucci et al. (2006) is also inconsistent with observations although the drop off is less extreme and can be brought into agreement with data by reducing the fraction of prompt SNe Ia from 50% to 30%. We conclude that the "prompt + delayed" description of the SN Ia rate must be only an approximation to a more smoothly declining DTD.

Zoom In Zoom Out Reset image size

Much of the stellar ejecta from CCSN, SN Ia, and AGB stars is not immediately returned to the cold ISM (see, e.g., Walch et al. 2015). Instead, it is injected into the warm ($T\gtrsim {10}^{4}$ K) ISM that does not form stars. Over time, gas in the warm ISM cools into the cold ISM and can get reincorporated into future generations of stars. While we have not included a warm ISM in the fiducial simulation, we explore its effects on chemical evolution in Figure 8.

Zoom In Zoom Out Reset image size

For the warm ISM simulation, we follow the scheme from Schönrich & Binney (2009a) and inject 99% of stellar yields from all sources (CCSN, SN Ia, and AGB stars) into the warm ISM with the remainder of the yields entering the cold ISM directly. In each time step, a fraction of the gas (dt/t_cool) cools out of the warm ISM and joins the cold ISM. We adopt a gas cooling timescale, t_cool, of 1.2 Gyr. The gas in the warm ISM is not ejected in outflows.

Including a warm ISM delays the return of yields from all sources, in contrast to changing the minimum time delay for SNe Ia (see Section 3.6), which only affects one yield source. The warm ISM simulation shows many signs of slower enrichment. It has much lower metallicity at early times, turns over at a lower [Fe/H], produces more high-α stars, and has broader main peaks to the MDF and [O/Fe]–DF than the fiducial simulation. Its equilibrium [O/Fe] and [Fe/H] are slightly higher than those of the fiducial simulation because the warm ISM prevents the ejection of a large fraction of the enrichment products, which is especially important for elements synthesized early, such as oxygen in CCSNe.

WAF derive analytic solutions for oxygen and iron evolution by adopting two simplifying assumptions beyond those used here: instantaneous recycling for material from AGB stars (but not for SN Ia enrichment), and supernova yields that are independent of metallicity. They assume an exponential or linear-exponential SFH and an exponential DTD for SNe Ia (or a sum of exponentials, which can be used to approximate the ${t}^{-1}$ power-law). Their analytic results capture many of the key features of Figure 3, including the dependence of equilibrium abundances and abundance ratios on η and SFE, the strong impact of SFE on the location of the knee in evolutionary tracks, the influence of the SFH (which is tied to the inflow rate), and the small impact of inflow metallicity. The analytic models also capture the dependence on the minimum delay time and shape of the DTD seen in Figures 5 and 7. These models do not incorporate a warm ISM, though they could be adapted to do so by introducing an exponential "DTD" for return of CCSN elements.

The adopted nucleosynthesis yields remain one of the largest sources of uncertainty in chemical evolution models (Timmes et al. 1995; Gibson 1997; Portinari et al. 1998; Romano et al. 2010). The yields of α-elements from CCSNe and iron-peak elements from SNe Ia are relatively robust, but the yields of elements that are sensitive to metallicity and elements produced in explosive burning phases or by neutron-capture are much more uncertain. Here we discuss our choice of yields and how they affect the chemical evolution model. We note that we use the term "yields" in the sense of solar masses of element X returned, as opposed to the classically defined sense of a mass fraction (e.g., Tinsley 1980).

4.1.1. CCSNe

4.1.2. SNe Ia

4.1.3. AGB Stars

In Figure 9, we compare the abundances from three simulations with different CCSN yields:

• 1.   
  Chieffi & Limongi (2004) and Limongi & Chieffi (2006) yields (the same as the fiducial simulation; hereafter CL04),
• 2.  
  Woosley & Weaver (1995, hereafter WW95) yields, and
• 3.  
  Limongi & Chieffi (2006, hereafter LC06) solar metallicity only yields.

We reduced the WW95 yields of the iron-peak elements (Cr, Mn, Fe, Co, Ni, Cu, and Zn) by a factor of two as recommended by Timmes et al. (1995) based on the normalization of [X/Fe]–[Fe/H] tracks. When we used the original WW95 yields, we found the expected ∼0.3 dex offset in [X/Fe] at low metallicity for non-iron-peak elements, though these yields were able to increase the final [Fe/H] by 0.1 dex to match the CL04 simulation. Timmes et al. (1995) note that this reduction is consistent with the uncertainty in the explosion energy that affects the mass cut and hence the yields of iron-peak elements. The WW95 iron yields have a metallicity dependence at low metallicity, which causes the rise in [O/Fe] as [Fe/H] decreases below −1. Our CL04 models, on the other hand, have metallicity independent iron yields by construction, so they produce nearly flat plateaus at low metallicity for elements such as O, Mg, Si, and S. Since the metallicity-dependence of the yields is a significant physical uncertainty, we also show the LC06 solar metallicity yields applied at all metallicities. The [X/Fe] trends for this simulation are often higher than the CL04 simulation, indicating that the solar metallicity CCSN yields are higher than the sub-solar metallicity CCSN yields (with the notable exceptions of Ca and Ti). Ni and Mn have similar SN Ia and solar metallicity CCSN yields, which produces their flat abundance trends. Elements with significant secondary production (e.g., N) show large discrepancies between the CL04 and LC06 trends. The Eu yields are metallicity independent by construction, so the differences between the CL04 and the solar metallicity CCSN yields show the slight metallicity dependence of the CL04 Fe yields.

Zoom In Zoom Out Reset image size

These simulations have a constant SFR but otherwise adopt the parameters of the fiducial simulation, including the SN Ia and AGB yields described above. Relative to the fiducial simulation, the constant SFR simulation has a lower equilibrium abundance and higher equilibrium [X/Fe] values (except Mn and Ni). Below we compare the predicted abundances of 20 elements for the two yield sets and discuss the agreement or disagreement with data from Reddy et al. (2003, 2006) and Ramírez et al. (2013) (oxygen only). Figure 10 shows the mass enrichment histories for these elements in the same format as Figure 2. Elements with metallicity-dependent yields have a steeper blue solid curve in Figure 10 than the red solid curve (iron mass). We discuss results element-by-element with reference to both figures. Our physical descriptions of the dominant enrichment sources and their metallicity dependence are based on our adopted yield models, which may or may not reflect reality.

Zoom In Zoom Out Reset image size

Carbon. Carbon enrichment, like oxygen, is dominated by massive stars. With our adopted yields, the contribution of AGB stars is sub-dominant but non-negligible above [Fe/H] ≈ −0.6 (see Figure 10). The [C/Fe] data follow an α-element trend except that the [C/Fe] turnover appears displaced toward somewhat higher [Fe/H]. The predicted [C/Fe] actually increases at [Fe/H] = −0.7, where carbon production by AGB stars peaks, which may explain this offset. At higher [Fe/H], AGB stars become net consumers of carbon. The CCSN C yields are mildly metallicity-dependent above $Z\sim 1/3\,{Z}_{\odot };$ however, they increase with Z for the CL04 yields but decrease for the WW95 yields. The CL04 simulation matches the observational data well for $[\mathrm{Fe}/{\rm{H}}]\gt -0.4$, but it predicts too low an abundance for the plateau (see the observed [C/Fe] trend from Bensby & Feltzing 2006, which the CL04 simulation agrees with fairly well). The WW95 simulation produces a similarly shaped trend to the data that is systematically too low by ∼0.3 dex. Interestingly, the LC06 solar metallicity yields reproduce the data very well.

Nitrogen. Most nitrogen production occurs in the CNO cycle, where carbon and oxygen get converted into nitrogen. Figure 10 shows that AGB stars dominate the production of nitrogen, especially for [Fe/H] < −1, so it is not surprising that the CL04 and WW95 simulations result in similarly shaped tracks in [N/Fe]–[Fe/H]. The [N/Fe] of the CL04 and WW95 simulations increases quickly from low metallicity to [Fe/H] ∼ −0.6 and −0.8, respectively, which reflects the secondary nature of nitrogen (i.e., its yield increases with the initial carbon and oxygen abundances). This metallicity dependent yield is also evident in the steepness of the blue dashed curve in Figure 10. It then declines as the nitrogen production from AGB stars is offset by losses to star formation and outflows while iron production from SNe Ia continues to increase (Figure 10). Both simulations underpredict [N/Fe] by ∼0.5 dex over the limited range of the data $-0.4\lesssim [\mathrm{Fe}/{\rm{H}}]\lesssim 0.1$.

Massive stars, (i.e., ones that will explode as CCSNe) could be an important source of nitrogen at zero and low metallicities. However, the uncertainty in the nitrogen yield of zero metallicity massive stars is believed to be at least an order of magnitude (Heger & Woosley 2010). These stars are born with low abundances of carbon and oxygen (relative to hydrogen), but rotational mixing could transport freshly synthesized carbon and oxygen from helium-burning layers to hydrogen-burning layers, making nitrogen production extremely sensitive to the strength of mixing. Nitrogen abundances are also difficult to measure due to the weakness of the Ni lines and their strong dependence on the measured effective temperature (Reddy et al. 2003). The [N/Fe]–[Fe/H] trend from the Reddy et al. (2003) data implies a super-solar [N/Fe] at solar metallicity, which suggests that either the Sun has a lower than typical [N/Fe] or, more likely, that the normalization of the Reddy et al. (2003) data is too high. Nonetheless, if observational estimates are approximately correct, then either the AGB or CCSN yields (or both) must be substantially in error.

Oxygen. Oxygen is produced almost exclusively by CCSNe (see Figure 10). Its yield has weak dependence on metallicity (see the LC06 curve) but a strong dependence on stellar mass, as can be seen from the effects of changing the slope or cutoff of the high mass end of the IMF in Figure 3(f). The parameters of the fiducial simulation were chosen to match the data in [O/Fe]–[Fe/H] and reach equilibrium at solar oxygen abundance and metallicity, so the CL04 simulation fits the data by construction. (While our simulation in Figures 9 and 10 assumes a constant SFR instead of exponentially declining infall, this has little impact on the shape of tracks as shown previously in Figure 3(b).) This model was not intended to fit the Ramírez et al. (2013) data above solar metallicity, as these stars were likely born interior to the solar annulus in the metal-rich inner disk and migrated to the solar neighborhood. The WW95 simulation fits the data well for [Fe/H] > −0.4, though it overpredicts [O/Fe] at lower metallicities ([O/Fe] declines from 0.85 to 0.6 from [Fe/H] = −2.0 to −1.0). Our model using the WW95 yields produces a 0.2–0.3 dex higher [O/Fe] plateau than Timmes et al. (1995) found with the same yields. While there are many differences between the two models, this discrepancy is almost certainly due to the fact that their model does not include yields for stars with M > 40 ${M}_{\odot }$ unlike our model (see Appendix B for more details). The differences between the WW95 and CL04 simulations at low metallicity are partially due to the metallicity dependence of the WW95 iron yields, whereas the CL04 iron yields are metallicity-independent by design since we chose a fixed mass cut. The range of plausible mass cuts and iron yields for CCSNe spans the range of uncertainty in the oxygen abundance measurements, which are challenging, even for the Sun (Caffau et al. 2011). In particular, abundances determined from the strong Oi triplet at 777 nm are sensitive to the adopted temperature scale and require corrections from the typical assumption of local thermodynamic equilibrium (LTE), adjustments of 0.1–0.5 dex depending on gravity and temperature (Ramírez et al. 2013).

Sodium. Nearly all sodium is synthesized in CCSNe. Its CCSN yield is strongly metallicity-dependent because it is an odd-Z element, hence the steepness of the blue dotted curve in Figure 10 and the rising trend of [Na/Fe] up to [Fe/H] < −0.4 for the CL04 simulation. At higher metallicities, [Na/Fe] decreases due to increased iron enrichment from SNe Ia. The CL04 simulation produces a dramatic rise in [Na/Fe] that is not seen in the observed [Na/Fe] trend, which is flat at [Na/Fe] ≈ 0.1 for $[\mathrm{Fe}/{\rm{H}}]\lt -0.4$ with a slow decline toward solar [Na/Fe] at solar [Fe/H]. However, low metallicity stars ([Fe/H] ≲ −1) can have sub-solar [Na/Fe] (e.g., Nissen & Schuster 1997, 2010), suggesting that the rising trend may qualitatively reflect observations at metallicities below that of the Reddy et al. (2003, 2006) data. The WW95 simulation initially declines from [Fe/H] = −2 to −1.6 due to its metallicity-dependent iron yields. It then rises from [Fe/H] = −1.6 to −0.6 but more gradually than the CL04 simulation due to a weaker metallicity dependence for the sodium yields. At [Fe/H] = −0.6, it turns over because of SN Ia iron enrichment. The normalization of the WW95 trend is a better match to the data than the CL04 simulation, and it does a reasonable job of fitting the data albeit with a sharper turnover. The LC06 trend always using the solar metallicity Na yield predicts an extremely high [Na/Fe] that is mostly off the top of the plot. There are some observational uncertainties in measuring sodium abundances, including non-LTE effects (see Smiljanic 2012 and references therein).

Magnesium. Like oxygen, magnesium enrichment is dominated by CCSNe, and the predicted yield is only mildly metallicity-dependent. The [Mg/Fe] data follow the typical α-element trend and are well reproduced by the CL04 simulation, though it slightly underpredicts [Mg/Fe], by ∼0.05–0.1 dex. The WW95 simulation also underpredicts [Mg/Fe], by ∼0.15–0.2 dex. Several authors (Timmes et al. 1995; Thomas et al. 1998; François et al. 2004) have found that magnesium is underproduced by the WW95 yields. François et al. 2004 noted that magnesium yields depend on the treatment of convection.

Aluminum. Aluminum is an odd-Z element produced almost entirely by CCSNe. The CL04 and WW95 simulations produce aluminum trends that are similar to those of sodium because the physics of production is similar. The CL04 simulation does a poor job of fitting the observational data for [Al/Fe], which instead resembles the α-element trend. The CL04 simulation overpredicts [Al/Fe] above [Fe/H] = −0.4 by ∼0.15 dex but dramatically underpredicts it at low metallicity. The weaker metallicity dependence of the WW95 yields relative to the CL04 yields results in a flatter [Al/Fe] trend that provides a better fit to the data. It underpredicts the normalization of the observations by ∼0.1 dex, but its shape is consistent with the data for [Fe/H] > −1. Though the normalization of the LC06 trend is too high, its shape matches the data well, adding further evidence that Al has a much weaker metallicity-dependence than the CL04 yields predict. Al abundance estimates can suffer non-LTE effects (see Andrievsky et al. 2008), especially at extremely low metallicities, which would affect any assessment of the low metallicity Al abundance trend and low metallicity CCSN yields.

Silicon. Silicon is an α-element made mostly in CCSNe like oxygen and magnesium but with a larger, though still small, contribution from SNe Ia. In CCSNe, it is produced in both hydrostatic and explosive burning phases (Woosley et al. 2002), and its yield is independent of metallicity. Both the CL04 and the WW95 simulations do an excellent job of reproducing the observed [Si/Fe] data. François et al. (2004) found good agreement with data when they adopted the solar metallicity WW95 silicon yields.

Sulfur. Sulfur is also an α-element, and its production is very similar to that of silicon with a slightly larger contribution from SNe Ia. Like silicon, sulfur is synthesized in hydrostatic and explosive burning phases (Woosley et al. 2002) of CCSNe, and it has a metallicity-independent yield. The CL04 and WW95 simulations agree with the mean of the [S/Fe] measurements, though these data span almost 0.5 dex in [S/Fe] and only 0.8 dex in [Fe/H].

Potassium. Potassium is an odd-Z element produced almost exclusively in CCSNe. Its yield is predicted to be metallicity-dependent, as can be seen from the super-linear slope of the blue dotted line in Figure 10, though its metallicity dependence is weaker than sodium or aluminum. The CL04 simulation significantly underpredicts the amount of potassium needed to achieve the solar value of [K/Fe] and the observed [K/Fe] data in Figure 9. The CL04 curve falls off the bottom of the plot, but the highest value it reaches is [K/Fe] = −0.65. The WW95 simulation also falls far below the data, though by 0.6 dex. The super-solar [K/Fe] values suggest that potassium may behave like an α-element at lower metallicities, but more data are needed to confirm such a trend.

Calcium. Calcium is an α-element made mainly in CCSNe, but with a significant SN Ia contribution. The CL04 simulation matches the level of the plateau in the [Ca/Fe] data, but its track turns over more sharply than the data and ends at a sub-solar [Ca/Fe] value. The CL04 CCSN yields decrease with metallicity, as shown by the LC06 track lying at lower [Ca/Fe]. The WW95 simulation produces higher [Ca/Fe] values at low metallicity, but it overlaps with the CL04 track for  $[\mathrm{Fe}/{\rm{H}}]\gt -0.8$. Mulchaey et al. (2014) argued that calcium-rich gap transients, a subclass of supernovae whose nebular spectra are dominated by calcium lines, are an important contributor of calcium in the intracluster medium, so it is possible that they might also be a meaningful source of calcium in galaxies.

Scandium. CCSNe are the dominant source of scandium. Like other odd-Z elements, the CL04 scandium yields have a very strong metallicity dependence that results in a rising trend with [Fe/H], which turns over at [Fe/H] ≈ −0.3. However, the observed scandium abundances more closely follow the typical α-element trend as in the metallicity-independent LC06 track. The WW95 simulation predicts an α-element trend for scandium, in stark contrast to the CL04 simulation. Both the CL04 and WW95 simulations underpredict the overall normalization of [Sc/Fe]. The offset for the WW95 simulation is ∼0.3–0.4 dex. A generic offset cannot be determined for the CL04 simulation because of its dramatically different shape compared to the data, though its peak is below the main trend by ∼0.15 dex.

Titanium. Titanium is both an α-element and a low iron-peak element. It is mostly produced in CCSNe, though SNe Ia contribute about 1/3 of the solar titanium abundance. Both the CL04 and the WW95 simulations produce α-element trends with similar normalizations. The data also follow an α-element trend. However, both simulations underpredict the titanium abundance data by ∼0.3–0.4 dex at all metallicities, and this underprediction is much more severe for the solar metallicity LC06 calculation. The underproduction of titanium is a generic problem for SN yields, and its cause is not well understood.

Vanadium. Vanadium is an odd-Z element produced predominantly in CCSNe but with a significant SN Ia contribution, so its predicted abundance trend behaves like a combination of odd-Z and Fe-peak elements. The CL04 vanadium yields have a weak metallicity dependence, and the CL04 simulation predicts a slightly rising [V/Fe] trend with [Fe/H] that peaks at [V/Fe] ≈ −0.4. However, the [V/Fe] data show a weak α-element trend that declines from [V/Fe] ∼0.2 to −0.1, which is roughly similar in shape to the LC06 curve but 0.3 dex higher in normalization. The WW95 simulation produces an α-element trend like the data, but it underpredicts the normalization by ∼0.4–0.5 dex.

Chromium. Chromium is an iron-peak element made equally in CCSNe and SNe Ia. The CL04 yields produce slightly super-solar [Cr/Fe], which declines to slightly sub-solar [Cr/Fe] due to the sub-solar [Cr/Fe] SN Ia yield ratio. The measured [Cr/Fe] data are centered on solar [Cr/Fe] with little scatter down to [Fe/H] = −1. The WW95 simulation declines by ∼0.2 dex from [Fe/H] = −2 to −1 and then rises slightly. It underpredicts the data by about 0.1 dex.

Manganese. SNe Ia and CCSNe contribute to the enrichment of manganese, an iron-peak element, at similar levels. The CL04 manganese yields increase with metallicity, resulting in a rising trend with [Fe/H]. The data show a unique observational trend of a plateau at [Mn/Fe] ≈ −0.4 for $[\mathrm{Fe}/{\rm{H}}]\lt -1$ that rises to solar [Mn/Fe] at solar metallicity. The CL04 simulation reproduces the data nearly perfectly, and the WW95 simulation reproduces the shape of the observed trend with a minor offset to lower [Mn/Fe]. The CCSN manganese yields are sensitive to the adopted mass cut, in a way that is not offset by similar changes to the iron yields. While the agreement between models and data in this panel is good, Bergemann & Gehren (2008) argue that the observed trend is caused by erroneously assuming local thermodynamic equilibrium, and that the trend of rising [Mn/Fe] toward higher metallicity goes away if non-LTE effects are taken into account.

Cobalt. Cobalt is an iron-peak element whose CCSN contribution is three times its SN Ia contribution. The CL04 cobalt yields are metallicity-dependent, which produces the rising [Co/Fe] trend with increasing [Fe/H]. In contrast to the CL04 simulation, however, the observed [Co/Fe] data follow an α-element trend, like that predicted by the solar-metallicity LC06 simulation but with a ∼0.2 dex offset. The WW95 simulation predicts a nearly identical trend to the CL04 simulation for $[\mathrm{Fe}/{\rm{H}}]\lt -0.8$, but it turns over at [Fe/H] = −0.6, whereas the CL04 simulation continues to rise until [Fe/H] = −0.3.

Nickel. Nickel is an iron-peak element that has slightly more contribution from SNe Ia than CCSNe. CCSN nickel yields are metallicity-dependent, so the CL04 simulation predicts a monotonically rising [Ni/Fe] trend. However, the observations show a flat trend at solar [Ni/Fe] with little scatter. The LC06 simulation produces a flat trend because the LC06 Ni/Fe ratio is similar to that of the SN Ia yields, though the normalization is 0.35 dex above the data. The WW95 prediction is similar to that of the CL04 simulation but offset toward higher [Ni/Fe] by ∼0.1–0.2 dex.

Copper. Copper is an iron-peak element that is exclusively produced in CCSNe. Both the CL04 and the WW95 copper yields have a strong metallicity dependence, and they predict rising [Cu/Fe] trends that turn over when SNe Ia begin to contribute significant amounts of iron. The [Cu/Fe] measurements roughly agree with this prediction, but the mean trend is offset to higher [Cu/Fe] by ∼0.1 dex. The data have significant scatter, especially at low metallicity. The WW95 simulation fits the low metallicity data better than the CL04 simulation because the former is offset toward lower [Fe/H].

Zinc. Zinc is an iron-peak element, whose yields are dominated by CCSNe. Both simulations produce an undulating trend that first declines, then rises due to a weak metallicity dependence of the zinc yields, and finally declines because of SN Ia iron enrichment. The observed [Zn/Fe] data show a weak α-element trend with large scatter below [Fe/H] = −1. The shape of the simulated trends may be consistent with the shape of the observed trend. However, the overall normalizations of the WW95 and CL04 predictions are too low by ∼0.25 dex and ∼0.5 dex, respectively, indicating that the CCSN yields are grossly underproducing zinc.

Europium. Europium is produced almost exclusively in the r-process. The astrophysical site of the r-process is unknown, but the two most likely candidates, CCSNe and neutron star mergers, are associated with massive stars. The [Eu/Fe] data follows an α-element trend, which suggests that europium yields are independent of metallicity. Neither the CL04 nor the WW95 yields included europium, so we adopted the Cescutti et al. (2006) Eu yields (their model 1), which are metallicity-independent by construction. The differences between the two simulations result from differences in the CCSN iron yields, and the small offset between the CL04 and the LC06 curves shows the slight metallicity dependence of the CL04 iron yields. The Cescutti et al. (2006) europium yields were determined empirically by fitting a chemical evolution model to observed [Eu/Fe] data from [Fe/H] = −3.4 to +0.15. Given the empirical nature of the europium yields, it is not surprising that the simulations and the data agree, but it is reassuring that we successfully reproduce a different observational [Eu/Fe] data set than the calibrating one with a different chemical evolution model.

For any individual element, the one-zone model predictions have freedom associated with the choice of SFE, outflow rate, and so forth, as discussed in Section 3. Once these parameters are chosen to produce a given [O/Fe] track, however, the predictions for other elements are sensitive mainly to yields. Thus, major discrepancies in Figure 9 probably indicate either incorrect data or incorrect yields. The most common forms of discrepancy are either an overall normalization offset or predictions of rising trends for elements with metallicity-dependent yields where observed trends are similar to those of α-elements. Adding more measurements at [Fe/H] < −1 would greatly increase the constraining power of the observational data.

Connecting the high-α sequence to the low-α sequence remains challenging for chemical evolution models. Motivated by early observational evidence of a gap in [Fe/O]–[O/H] data from Gratton et al. (2000); Chiappini et al. (2001) found that a track could produce a gap in the abundance trend if the disk formed in two major infall episodes separated by a cessation of star formation. The gap in the SFH is caused by a slow infall rate for the second episode coupled with a threshold gas density for star formation (Kennicutt 1989). The effect of this threshold is magnified in the Chiappini et al. (2001) model because it increases from 4 to 7 ${M}_{\odot }$ pc⁻² between the first and second infall episodes. Figure 11(a) shows a pair of simulations of the two-infall model for hiatuses in star formation of 200 Myr (blue points) and 1 Gyr (red points). The points show stars randomly chosen from the simulations with Gaussian noise of σ = 0.05 in [Fe/H] and σ = 0.02 in [O/Fe] added. For the first Gyr, we set the inflow timescale to be 1 Gyr and set the SFE to be $6\times {10}^{-10}\,{\mathrm{Gyr}}^{-1}$. Next, we completely suppressed star formation starting at t = 1.0 Gyr for 200 Myr or 1 Gyr. Then, we set the inflow timescale to be 6 Gyr and reduced the SFE by a factor of two (to $3\times {10}^{-10}$ ${\mathrm{Gyr}}^{-1}$) with a maximum inflow (${\tau }_{\max }$) at 1 Gyr ($\propto {e}^{-(t-{\tau }_{\max })/{\tau }_{1}}$).

Zoom In Zoom Out Reset image size

The inflow timescale for the first epoch and the ratio of the SFE between the epochs (a factor of 2) were adopted from the original formulation of the two-infall model in Chiappini et al. (1997), but absolute values of SFE are not easily comparable because they adopted a different star formation law. However, the exact values of these parameters are much less significant than the choice of the length of the hiatus in star formation. Our choices of hiatuses of 200 Myr and 1 Gyr are intended to demonstrate the sensitivity of the abundances to this parameter.

The two-infall model produces a gap in [O/Fe]–[Fe/H] because of the hiatus in star formation. During the hiatus, SN Ia iron enrichment decreases [O/Fe], but its effect of increasing [Fe/H] is offset almost perfectly by dilution from inflow. The magnitude of the gap depends sensitively on the length of the hiatus in star formation. For a 200 Myr hiatus, [O/Fe] drops by 0.07 dex, and the subsequent evolutionary track is roughly a continuation of the pre-hiatus track but with a shallower slope. For a 1 Gyr hiatus, [O/Fe] drops 0.23 dex to [O/Fe] ≈ 0.0. However, the sudden return of star formation boosts the rate of CCSNe relative to SNe Ia, which quickly drives [O/Fe] back up to $\sim +0.1$, after which the model evolves to its solar abundance equilibrium. In both models, the pause in star formation leads to a separate peak in the MDF at [Fe/H] ≈ −0.3. The large gap, inverted-U shape of the low-α sequence, and the bimodal low-α sequence of the 1 Gyr hiatus simulation are not seen in the observational data, implying that the length of the hiatus must fall in a narrow range around 200 Myr to produce the right size gap and the morphology of the low-α sequence. However, we have not found any combination of parameters that reproduces the morphology of the observed stellar distribution, with two distinct sequences showing high and low [α/Fe] over a substantial range in [Fe/H].

Schönrich & Binney (2009a, 2009b) incorporated the effects of radial mixing of gas and stars into a chemical evolution model. They included the effects of changing the epicyclic amplitude ("blurring") and guiding center radius (radial migration or "churning"). Their model produces the high-α and low-α sequences from a superposition of stars born at a range of radii, without a hiatus in star formation. The high-α sequence is formed by overlapping tracks from different radii. Once SNe Ia become significant contributors of iron, the tracks move quickly across the valley in [α/Fe]–[Fe/H]. Like all models with continuous SFHs, their model produces some but not many intermediate [α/Fe] stars. The low-α sequence is a superposition of the stars near the equilibrium abundances of a range of radii, not an evolutionary track. Since individual tracks rapidly asymptote to their equilibrium abundance, the low-α sequence contains many more stars than the valley at intermediate [α/Fe].

Figure 11(b) shows a superposition of simulations that was designed to reproduce the high-α and low-α sequences by varying the outflow mass-loading parameter and inflow timescale to mimic the enrichment histories of several different galactocentric radii, whose stars were then mixed together. Increasing η decreases the equilibrium [Fe/H] (see Figure 3(d)) and results in a range of equilibrium metallicities that form a ridge line at low-α. Increasing the inflow timescale maintains the same trajectory of a simulation but moves its equilibrium abundance higher up the track to higher [O/Fe] and lower [Fe/H], which produces the negative slope of the low-α sequence. To reproduce the slope, we adopted a constant SFR for the highest-η simulation that effectively acts as an infinite inflow timescale. As the inflow timescale and η increase along our model sequence, the final stellar mass of the models decreases (by a factor of three from the lowest-η to the highest-η simulations) because the simulations accrete less gas and retain that gas less efficiently. Our choices of outflow mass-loading parameter and inflow timescale are intended to be physically plausible values that qualitatively reproduce a two-sequence abundance trend.

The points in Figure 11(b) are stars randomly drawn from the simulations with Gaussian noise of $\sigma =0.05$ in [Fe/H] and $\sigma =0.02$ in [α/Fe] added to show the 2D distribution of stars in [O/Fe]–[Fe/H]. The composite [O/Fe]–DF (black line) shows a single peak at low-α, in contrast to the two-infall model with a 1 Gyr hiatus. To account very roughly for the fact that stars born in situ are more common than migrated stars, we down-selected stars from the simulations by a factor of 0.2, 1. (fiducial), 0.2, and 0.1 (from lowest to highest η). This suite of simulations is intended to illustrate how the superposition scenario works, and a more complete multi-zone model with an accurate treatment of stellar migration and radial gas flows (such as Schönrich & Binney 2009a, 2009b) is needed to make more quantitative comparisons with data.

While these parameter variations were motivated by the results shown in Figures 3(a) and (d), they reflect differences in SFH as a function of galactocentric radius. An increasing inflow timescale with radius is a generic aspect of inside-out galaxy formation (Larson 1976). The outflow mass-loading parameter also might increase with radius because the outer disk may be less effective at retaining gas and metals due to a weaker vertical potential and a lower gas density to counteract energy injection by SNe. Furthermore, gas and metals flow inward through the disk (Stark 1984), enriching the inner disk at the expense of the outer disk, which mimics the effect of the outflow mass-loading parameter increasing with radius. While the superposition model of Figure 11(b) does have high-α and low-α stars at the same [Fe/H], the bimodality of the distribution is arguably weaker than observed (though one must be cautious of observational selections that may amplify this bimodality). The low-α ridge line could be extended to lower [Fe/H] by adding models with still higher η, representing migration from the far outer disk.

To summarize, a single simple chemical evolution model with constant parameters cannot reproduce both the high-α and low-α sequences. The two-infall model produces a gap between the high-α and low-α sequences due to a gap in the SFH. The length of this hiatus in star formation must be narrowly confined to about 200 Myr to avoid overshooting the low-α sequence. Alternatively, the high-α and low-α sequences can be formed from a mix of stellar populations, born at different galactocentric radii with distinct enrichment histories, that migrated to the solar neighborhood.

WAF show that one can produce a bimodal [α/Fe]-[Fe/H] distribution in a one-zone model by sharply increasing η after the model has evolved to an initial equilibrium, so that the stellar population then evolves backward to lower [Fe/H] at roughly constant [α/Fe]. We have not examined such a scenario here, but evolution of η, radial mixing, and unsteady star formation may all play a role in shaping the observed distribution of stars in [α/Fe]–[Fe/H] space.

The left panels of Figure 12 shows the eigenvector components for the first two PCs of the multi-element abundance space for the fiducial simulation. The points represent the contribution of the elemental abundances relative to Fe for each PC. A large eigenvector component (positive or negative) indicates that an elemental abundance contributes heavily to a PC, whereas an eigenvector component close to 0 implies that the PC is nearly independent of that elemental abundance. Two elements have correlated contributions to a PC if their eigenvector components have the same sign and anticorrelated contributions if their eigenvector components have opposite signs. The eigenvector components are normalized such that they quadratically sum to unity. In Figure 12, we show PC1 and PC2 for the full sample (black) and for a population of stars with an equal number sampled in 0.1 dex bins of [Fe/H] from −1.0 to +0.1 (gray). The equal metallicity sampling was designed to roughly represent selection biases in existing multi-element surveys such as Reddy et al. (2003, 2006) and Bensby et al. (2003, 2005, 2014).

Zoom In Zoom Out Reset image size

PC1 of the full sample is dominated by a correlation among the abundances of α-elements (C, O, Mg, Si, S, Ca, and Ti) and N. These abundances show a slight anticorrelation with the abundances of Mn and Ni, two Fe-peak elements with metallicity-dependent CCSN yields. PC1 is mostly independent of the abundances of the odd-Z elements (Na, Al, and K), and the Fe-peak elements, V, Cr, and Co. Thus, PC1 mainly captures the relative contribution of CCSNe and SNe Ia, the primary factor in multi-element abundance evolution.

PC2 of the full sample is characterized by a strong correlation among the abundances of N, the odd-Z elements, and the Fe-peak elements V, Mn, Co, and Ni. PC2 is independent of α-element abundances. The main contributing elemental abundances to PC2 all have rising trends at low metallicity because of metallicity-dependent CCSN or AGB star (N) yields.

When we apply PCAA to the equal metallicity sampling, however, PC1 and PC2 are swapped. PC1 of the equal metallicity sampling shows an anticorrelation between elements with metallicity-dependent yields (odd-Z elements, Mn, and Ni) and α-elements. Instead of tracking the decline from the knee due to SNe Ia, it tracks the progression along the high-α plateau at low [Fe/H], like PC2 of the full sample. PC2 of the equally metallicity sampling is a correlation among all non-Fe-peak elements, which shows the dilution of the abundances of these elements by Fe production from SNe Ia. This PC is similar to PC1 of the full sample, except that PC1 of the equally metallicity sampling has already accounted for the odd–even effect.

The right panel of Figure 12 shows the distribution in PC1–PC2 space of randomly selected stars (analogous to the full sample) with Gaussian noise of σ = 0.02 in both PC1 and PC2 added to display overlapping points. Each cluster of points is a stellar generation, and these clusters overlap as the simulation approaches the equilibrium abundances. The abundances evolve in PC1–PC2 space from low to high PC2 values in the first few hundred Myr and then move toward lower PC1 values over the remainder of the simulation. The initial increase in PC2 values at early times is driven by the rapid increase in the abundances of elements with metallicity-dependent yields. The PC1 values decrease as SN Ia Fe enrichment (starting around [Fe/H] = −0.8) decreases the relative abundances of elements with metallicity-independent CCSN yields (C and α-elements) and N, whose AGB star yields are less metallicity-dependent than at lower metallicities.

The lack of contribution, positive or negative, to PC1 from the abundances of the odd-Z elements Na, Al, K, and V is a consequence of the offset between the knee in the abundances of α-elements (plus C and N) at [Fe/H] ≈ −1 and the turnover in the abundances of odd-Z elements at [Fe/H] ≈ −0.4. The abundances of odd-Z elements continue to increase after the knee in α-element abundances. The abundances of both α- and odd-Z elements decline after the turnover in odd-Z element abundances. These two trends cancel out because the PCs are linear by construction. A nonlinear dimensionality reduction technique would be better suited to this case and possibly in general given the nonlinearity of nucleosynthesis processes, but we defer any such explorations to future work.

We applied PCAA to each of the parameter variation simulations in Figures 3, 5, 7, and 8. The PC1s of the full sample in these simulations are generally similar to PC1 of the fiducial simulation and all have significant contributions from α-elements. The contribution to PC1 from elements with metallicity-dependent CCSN yields can be correlated, uncorrelated, or anti-correlated with the α-elements depending on

• 1.  
  the offset in metallicity of the α-element knee and the turnover of the metallicity-dependent elemental abundance trends, and
• 2.  
  the magnitude of the decline in the metallicity-dependent elemental abundance trends after the turnover.

A smaller offset between the knee and the turnover produces a correlation between α- and metallicity-dependent elements, such as in the high SFE (higher metallicity knee) and the high outflow mass-loading parameter (lower metallicity turnover) simulations. These two simulations also feature a steep decline in the metallicity-dependent elemental abundances above the turnover, which contributes to the correlation. On the other hand, simulations that have nearly monotonically rising metallicity-dependent elemental abundance trends—like the low SFE, population synthesis SN Ia DTD, and warm ISM simulations—show a strong anti-correlation between α- and metallicity-dependent elements.

The difference between full sampling and equal metallicity sampling is a caution that PCs depend on not just the enrichment physics but the sampling of stars. This is unsurprising, as it is capturing the directions of greatest variance, and the sampling will affect which variations have the biggest quantitative effect. When comparing models to data one therefore needs to match the sample properties. However, with fixed sampling the structure of the PCs does depend on the chemical evolution physics. For subsequent applications of PCAA, we will use the full sample unless specified otherwise.

6.2.1. CL04 versus WW95 CCSN Yields

Figure 13 shows the effect of varying the CCSN yields on the PC1 and PC2 eigenvector components. The black and gray lines represent the CL04 yield (same as Figure 12) and the WW95 yield simulations, respectively. Both simulations were run with the fiducial parameters. PC1 of the WW95 simulation shows a correlation between the abundances of α- and odd-Z elements because of the weaker odd–even effect in the yields. The WW95 PC1 describes more of the variation than the CL04 PC1 (89.6% versus 68.4%). Like the CL04 PC2, the WW95 PC2 captures the correlated abundances of N, Na, Al, Co, and Ni, though it is nearly independent of the abundances of K, V, and Mn and shows a weak anticorrelation with the heavier α-element abundances. This comparison demonstrates that PC structure is a potential probe of SN yield physics provided that the CCSN contribution can be separated from other effects.

Zoom In Zoom Out Reset image size

6.2.2. IMF Variations

Figure 14 shows PC1, PC2, and PC1–PC2 space for PCAA of a 1:1 mixture of the Kroupa and flat (α = 2.1) IMF simulations. We applied PCAA to this mixture of simulations as a conceptual model of the mixture of two stellar populations formed with different IMFs that might result from the accretion of satellites or radial mixing in the disk. PC1 and PC2 mostly resemble PC1 and PC2 of the fiducial simulation with some key differences, such as a correlation among odd-Z and α-elements in PC1. More importantly, the two sub-populations are cleanly separated in PC1–PC2 space. While their tracks in [O/Fe]–[Fe/H] are also distinct (see Figure 3(f)), PCAA coherently adds the differences in all elements to amplify small offsets between populations, which could be used to identify more subtle IMF variations between populations with otherwise similar formation histories.

Zoom In Zoom Out Reset image size

In Figure 15, we show PC1–PC2 space for the superposition scenario discussed in Section 5.2. PC1 and PC2 (and PC3) of the superposition scenario are nearly identical to those of the fiducial simulation (left panels of Figure 12) for either the full sample or equal metallicity sampling. The superposition scenario differs not in the structure of the PCs but in the distribution of stars in the space. The four simulations follow the same path in PC1–PC2 space at early times (high PC1 and low PC2) but separate at late times (low PC1 and high PC2). This separation is cleaner than in [O/Fe]–[Fe/H] (Figure 11(b)) and is particularly noticeable in PC2 because of its large contributions from elements with metallicity-dependent yields. Consequently, the fractions of the variation attributed to PC1 and PC2 are slightly smaller and larger, respectively, than those of the fiducial simulation. This scenario shows that PCAA can separate stellar populations with distinct histories, using empirical correlations to take advantage of all the information in the data, in contrast to focusing on a single element or group of elements to define populations.

Zoom In Zoom Out Reset image size

In Figure 16, we show PC1, PC2, PC3, and PC4 for the observed abundances of Reddy et al. (2003, 2006) and Ramírez et al. (2013) shown in Figure 9 (which we will refer to as the Reddy–Ramírez data). We adopted the superposition scenario as the best model comparison for the Reddy–Ramírez data because it more closely reproduces the observed scatter in abundances than one-zone simulations. Additionally, the gray solid line shows a down-selection of stars from the superposition scenario to match the Reddy–Ramírez iron MDF. We note that the abundances of N, S, and K were not included in the analyses discussed in this section for either the superposition scenario or the observed abundances due to limited data.

Zoom In Zoom Out Reset image size

PC1 of the Reddy–Ramírez data broadly matches the α-element abundance nature of PC1 of the metallicity-matched and full samples of the superposition scenario. However, the abundances of the odd-Z elements Na, Al, and V and the Fe-peak elements Co and Ni are correlated with the abundances of α-elements. This is evidence that the metallicity dependence of the CL04 yields is too strong (the metallicity-dependence for these elements is weaker for the WW95 yields) or else that delayed gas mixing can significantly dampen its effects. The abundance of Al contributes more strongly to PC1 of the data than the abundance of Na, whereas the abundances of these elements generically contribute to PC1 and PC2 of the simulations in lockstep for both the CL04 and WW95 yields. This suggests that the metallicity dependence of Na and Al CCSN yields is more distinct than the CL04 and WW95 yields predict, though it is also possible that non-LTE effects are biasing the observed abundances of Na and Al. Both the CL04 and WW95 yields overpredict the metallicity dependence of Co and Ni.

PC2 of the data is dominated by the abundance of C, in contrast to PC2 of the metallicity-matched and full samples of the superposition scenario, which are driven by the abundances of elements with metallicity-dependent CCSN yields. The overwhelming importance of the C abundance to PC2 suggests that PC2 could reflect variation due to AGB star nucleosynthesis, either in the standard sense of a well-mixed gas reservoir or in binary mass transfer systems, or observational scatter in the C abundances. We note that PC2 of the data accounts for a much lower fraction of the variation (9% versus 26%) than PC2 of the metallicity-matched superposition scenario, so any interpretation of its potential physical origins should be treated with caution.

PC3 describes a comparable fraction of the variation as PC2 (7%). It is dominated by Mn and the other Fe-peak elements, tentatively hinting at variations in the SN yields of these elements. Alternatively, and less interestingly, PC3 might result from variance caused by random observational errors in the Fe abundance, shifting the observed Fe coherently relative to these other elements that are tied closely to the true Fe abundance.

PC4 also accounts for a non-negligible fraction of the variation (5%). However, it shows Co anticorrelated with O and V, whose physical interpretation is unclear and is likely not warranted. We conjecture that this PC reflects mainly observational noise and limited sample size.

Figure 17 shows PC1–PC2 space for the Reddy–Ramírez data. Metallicity generally increases from right to left. The data exhibit a bimodal distribution in PC1, reflecting the bimodality between the high-α and low-α sequences. However, the kinematic selection of thick disk stars in this sample has enhanced this bimodality to create two distinct populations. In contrast to the superposition scenario (see Figure 15), PC2 shows no trend with PC1 and substantial scatter. However, the nature of PC2 is very different in these two cases.

Zoom In Zoom Out Reset image size

Figure 18 shows the cumulative fraction of variation attributed to the first N PCs for the Reddy–Ramírez data (black solid line), the superposition scenario sampled to match the Reddy–Ramírez MDF (gray solid line), and the full superposition scenario sample (gray dashed line). The first two PCs of the superposition scenario explain almost 99% of the variation, while eleven PCs are required to reach the same threshold for the data. The much larger number of PCs required to explain the data could be a consequence of more complex physics or more stochastic enrichment affecting the data, or it could be a consequence of observational errors producing variance, or both. To fully explore the contribution of observational errors, one should create synthetic spectra from the model distribution of stars, add noise, and analyze them like one does the data. This large undertaking is beyond the scope of this paper. An alternative to fraction of variance is to ask what fraction of stars are well fit in a ${\chi }^{2}$ sense by 1 PC, 2 PCs, ... N PCs, or to plot the distribution of ${\chi }^{2}$ values for stars after fitting with 1 PC, 2 PCs, etc. These quantifications should be more robust to observational noise, at least if the noise itself is accurately determined.

Zoom In Zoom Out Reset image size

To summarize the PCAA data–model comparisons:

• 1.  
  The models and data both have the abundance of CCSN elements relative to Fe-peak elements as the strongest PC.
• 2.  
  The models with our fiducial CCSN yields (CL04) predict a strong signature of metallicity-dependent yields for Na and Al (also for N and K but these elements have limited data). This trend is weaker in the data, particularly for Al. It would be very interesting to have N and K observations as well.
• 3.   
  The models predict that C behaves similarly to other α-elements, but the data show additional variation in C essentially on its own as PC2. This could result from enrichment physics, binary physics, or observational errors.
• 4.  
  The data have a correlated contribution of Fe-peak elements relative to Fe, though it is not clear whether this is physically significant or an observational effect.
• 5.  
  The data require many more PCs to explain 99% of the variance. Observational noise certainly contributes, but more detailed study is needed to decide whether it is the dominant factor or whether the larger number of PCs reflects more complex enrichment physics.

We start by initializing the model with the initial gas mass with primordial mass fractions. At each time step, we update the mass of each isotope, i, as follows,

$$\begin{eqnarray}{M}_{i}({t}_{j}) & = & {M}_{i}({t}_{j-1})+{M}_{i,\mathrm{inflow}}(t)-{M}_{i,\mathrm{SF}}(t)-{M}_{i,\mathrm{outflow}}(t)\\ & & +{\mathrm{CCSN}}_{i}(t)+{\mathrm{AGB}}_{i}(t)+\mathrm{SN}\,{\mathrm{Ia}}_{i}(t),\end{eqnarray} \tag{ 9 }$$

where ${M}_{i,\mathrm{inflow}}$ is the inflowing mass, ${M}_{i,\mathrm{SF}}$ is the mass removed by star formation, ${M}_{i,\mathrm{outflow}}$ is the outflowing mass, and CCSN_i, AGB_i, and SN Ia_i are the absolute yields from each of these sources. We describe each of these terms in more detail below.

The inflow mass of a given isotope is determined by multiplying the inflow rate function (${\dot{M}}_{\mathrm{inflow}}(t)$), the mass fraction of the isotope in the inflowing gas (${f}_{i}^{\mathrm{inflow}}$), and the length of the time step ${\rm{\Delta }}t$:

$$\begin{eqnarray}&&{M}_{i,\mathrm{inflow}}(t)={\dot{M}}_{\mathrm{inflow}}(t)\cdot {f}_{i}^{\mathrm{inflow}}\cdot {\rm{\Delta }}t.\end{eqnarray} \tag{ 10 }$$

In the fiducial model the inflow mass fractions are primordial. The inflow rate can be an exponential (Equation (1)), a double exponential (Equation (2)), a linear–exponential product (Equation (3)), or set to produce a constant SFR by maintaining a constant total gas reservoir mass by exactly offsetting the other terms in Equation (9).

One of the unique features of our model is that it does not rely on population averages but stochastically forms every individual star. In each time step, the number of stars born in each mass bin (of width 0.1 ${M}_{\odot }$ from 0.1 to 8.0 ${M}_{\odot }$ and 1.0 ${M}_{\odot }$ from 8.0 to 100 ${M}_{\odot }$) is determined by the following steps.

• 1.  
  Calculate the expectation value of the SFR (SFR^EV) from the Kennicutt–Schmidt law,

  $$\begin{eqnarray}&&{\mathrm{SFR}}^{\mathrm{EV}}(t)=\displaystyle \frac{1}{{t}_{\mathrm{gas}}}{Area}{\left(\displaystyle \frac{{M}_{\mathrm{gas}}}{{Area}}\right)}^{N},\end{eqnarray} \tag{ 11 }$$

  where N is the Kennicutt–Schmidt law index. In the fiducial model, we adopt N = 1, so the area terms cancel and Equation (11) reduces to Equation (4).
• 2.  
  Multiply SFR^EV by the length of the time step to get the expectation value of the mass of newly forming stars (SF^EV),

  $$\begin{eqnarray}&&{{SF}}^{\mathrm{EV}}(t)={\mathrm{SFR}}^{\mathrm{EV}}(t)\cdot {\rm{\Delta }}t.\end{eqnarray} \tag{ 12 }$$
• 3.  
  Distribute SF^EV into mass bins according to the IMF, so the expectation value of the mass of newly formed stars in bin m_j (${\mathrm{SF}}_{{m}_{j}}^{\mathrm{EV}}$) is

  $$\begin{eqnarray}&&{\mathrm{SF}}_{{m}_{j}}^{\mathrm{EV}}(t)={\mathrm{SF}}^{\mathrm{EV}}(t)\left(a{\int }_{{m}_{\mathrm{bin}\mathrm{low}}}^{{m}_{\mathrm{bin}\mathrm{up}}}{m}^{-\alpha +1}{dm}\right),\end{eqnarray} \tag{ 13 }$$

  where ${m}_{\mathrm{bin}\mathrm{low}}$ and ${m}_{\mathrm{bin}\mathrm{up}}$ are the lower and upper mass limits of mass bin m_j, a is the appropriate normalization constant for the chosen IMF, and α is the IMF power law slope over the mass bin.
• 4.  
  Convert ${\mathrm{SF}}_{{m}_{j}}^{\mathrm{EV}}$ into the expected number of stars in mass bin m_j (${N}_{\mathrm{SF},{m}_{j}}^{\mathrm{EV}}$) by dividing by the average mass of the bin:

  $$\begin{eqnarray}&&{N}_{\mathrm{SF},{m}_{j}}^{\mathrm{EV}}(t)={\mathrm{SF}}_{{m}_{j}}^{\mathrm{EV}}(t)\displaystyle \frac{{\int }_{{m}_{\mathrm{bin}\mathrm{low}}}^{{m}_{\mathrm{bin}\mathrm{up}}}{m}^{-\alpha }{dm}}{{\int }_{{m}_{\mathrm{bin}\mathrm{low}}}^{{m}_{\mathrm{bin}\mathrm{up}}}{m}^{-\alpha +1}{dm}}.\end{eqnarray} \tag{ 14 }$$
• 5.  
  Stochastically sample the expected number of stars in mass bin m_j according to a Poisson distribution to get the actual number of stars formed in each mass bin (${N}_{\mathrm{SF},{m}_{j}}$):

  $$\begin{eqnarray}&&{N}_{\mathrm{SF},{m}_{j}}(t)=\mathrm{Poisson}({N}_{\mathrm{SF},{m}_{j}}^{\mathrm{EV}}(t)).\end{eqnarray} \tag{ 15 }$$

  For the fiducial model, the effect of random sampling of the IMF is small because of its large mass, so the maximum difference between the expected value of the SFR and the actual SFR was only 0.06%.
• 6.  
  Compute the actual total mass in newly formed stars by multiplying the number of stars formed in each mass bin by the average mass in the bin and then summing over all mass bins:

  $$\begin{eqnarray}&&\mathrm{SF}(t)=\displaystyle \sum {\mathrm{SF}}_{{m}_{j}}(t)\\ &&\quad =\displaystyle \sum {N}_{\mathrm{SF},{m}_{j}}(t)\displaystyle \frac{{\displaystyle \int }_{{m}_{\mathrm{bin}\mathrm{low}}}^{{m}_{\mathrm{bin}\mathrm{up}}}{m}^{-\alpha +1}{dm}}{{\displaystyle \int }_{{m}_{\mathrm{bin}\mathrm{low}}}^{{m}_{\mathrm{bin}\mathrm{up}}}{m}^{-\alpha }{dm}}.\end{eqnarray} \tag{ 16 }$$
• 7.  
  Finally, multiply the total mass in newly formed stars by the mass fraction of isotope i in the ISM at time t (f_i(t)) to get the mass of isotope i that goes into newly formed stars,

  $$\begin{eqnarray}&&{\mathrm{SF}}_{i}(t)=\mathrm{SF}(t)\cdot {f}_{i}(t).\end{eqnarray} \tag{ 17 }$$

In the fiducial model, the outflow rate is a constant factor times the SFR. The composition of the outflowing gas is the same as the ISM composition, so the mass of isotope i ejected from the galaxy in outflows is simply:

$$\begin{eqnarray}&&{\mathrm{outflow}}_{i}(t)=\eta \cdot {\mathrm{SF}}_{i}(t),\end{eqnarray} \tag{ 18 }$$

where η is the outflow mass-loading parameter.

For each time step, we determine the mass range of stars from each previous time step that will evolve in the current time step by inverting the stellar lifetime function. If only some of the stars in a mass bin will evolve because the age range of the mass bin straddles the edge of the time step, then the model calculates the fraction of stars in that mass bin that will evolve. This step is particularly helpful in smoothing out the gas return from very long-lived stars. We then look up the appropriate CCSN or AGB star yields for each mass bin and the metallicity of the ISM at birth to compute the yields. We also add the newly formed remnant mass to the existing remnant population.

To calculate the mass return from SNe Ia, we use the SN Ia DTD to compute the expectation value of the number of SNe Ia that will explode in a given time step and then stochastically sample from a Poisson distribution to determine the actual number of SNe Ia. The calculation of the expected number of SNe Ia varies among the four different DTDs that we considered (exponential, power law, population synthesis, and prompt + delayed), so we will describe the unique aspects of the procedures for each DTD separately.

A.6.1. Exponential SN Ia DTD

A.6.2. Power-law and Population Synthesis SN Ia DTDs

For the power-law and population synthesis DTDs, we first compute the SN Ia rate per solar mass of star formation, given by Equation (7) for the power-law DTD and described in detail in Section 3.1 of Greggio (2005) for the population synthesis DTD. Next, we normalize the SN Ia rate to match observational constraints of the total number of SNe Ia from 40 Myr to 10 Gyr of an individual stellar population (see Section 3.6 for details). For previous stellar generations, we calculate the product of the SN Ia rate and the stellar mass formed in that time step (SF), which gives us the expected number of SNe Ia (${N}_{\mathrm{SN}\,\mathrm{Ia}}^{\mathrm{EV}}({t}_{1})$) from stars born in that time step:

$$\begin{eqnarray}&&{N}_{\mathrm{SN}\,\mathrm{Ia}}^{\mathrm{EV}}({t}_{1})={R}_{\mathrm{SN}\mathrm{Ia}}({t}_{1})\cdot \mathrm{SF}({t}_{1}).\end{eqnarray} \tag{ 21 }$$

Then, we sum over the contributions of the previous time steps older than the minimum delay time to get the total expected number of SNe Ia (${N}_{\mathrm{SN}\,\mathrm{Ia}}^{\mathrm{EV}}$):

$$\begin{eqnarray}&&{N}_{\mathrm{SN}\,\mathrm{Ia}}^{\mathrm{EV}}=\displaystyle \sum _{{t}_{1}=0}^{{t}_{1}=t-{t}_{\min }}{N}_{\mathrm{SN}\,\mathrm{Ia}}^{\mathrm{EV}}.\end{eqnarray} \tag{ 22 }$$

A.6.3. Prompt + Delayed DTD

The expected number of SNe Ia for the prompt + delayed DTD is simply the SN Ia rate as expressed in Equation (8) multiplied by the length of the time step (for all time steps after the minimum delay time):

$$\begin{eqnarray}&&{N}_{\mathrm{SN}\,\mathrm{Ia}}^{\mathrm{EV}}(t)={R}_{\mathrm{SN}\mathrm{Ia}}(t)\cdot {\rm{\Delta }}t.\end{eqnarray} \tag{ 23 }$$

A.6.4. Stochastic Realization of SNe Ia

Once we have the expectation value of the number of SNe Ia to explode, we compute the actual number of SNe Ia by drawing from a Poisson distribution:

$$\begin{eqnarray}&&{N}_{\mathrm{SN}\mathrm{Ia}}(t)=\mathrm{Poisson}({N}_{\mathrm{SN}\,\mathrm{Ia}}^{\mathrm{EV}}(t)).\end{eqnarray} \tag{ 24 }$$

The mass of isotope i returned to the ISM from SNe Ia is simply the number of SNe Ia multiplied by the absolute SN Ia yield of isotope i (${Y}_{\mathrm{SN}\mathrm{Ia},{\rm{i}}}$) in solar masses,

$$\begin{eqnarray}&&\mathrm{SN}\,{\mathrm{Ia}}_{i}(t)={N}_{\mathrm{SN}\mathrm{Ia}}(t)\ast {Y}_{\mathrm{SN}\mathrm{Ia},i}.\end{eqnarray} \tag{ 25 }$$

Finally, we subtract the total mass returned from SNe Ia from the total mass in remnants (for all DTDs) and from reservoir of WDs that will explode as SNe Ia (i.e., ${M}_{\mathrm{WD}}^{\mathrm{SN}\,\mathrm{Ia}}$) (for the exponential DTD only).