12. Chemical evolution in galaxies¶
So far, we have mainly focused on the dynamics of galactic systems, but galaxies are also the stage on which stars and planets form, evolve, and die. The formation and evolution of galaxies and the stars within them are intimately connected and in detail it is impossible to study the evolution of either one of them without considering the other. In this chapter, we start looking at the astrophysical processes important for galaxy formation and evolution in more detail, by considering how the abundances of different chemical elements in interstellar gas and stars evolves over time.
12.1. Chemical enrichment processes in galaxies¶
On a basic level, galaxies consist of dark matter, stars, and gas. After the formation of the dark-matter halo (discussed in more detail in Chapter 18), dark matter plays a largely passive role as the gravitational stage upon which stars and gas perform an elaborate dance. Stars in the Universe by and large form in the disks of large galaxies like the Milky Way, from cold gas that forms a thin layer embedded within the thicker galactic disk (the interstellar medium; ISM). When stars form from gas that has cooled to form molecular clouds, they inherit the abundances of the ISM at the time of their formation. As stars age and massive stars die, some or all of their material is returned to the ISM as ejecta, (partially) enriched through nucleosynthesis during a star’s regular life and (sometimes) violent, explosive death. This ‘enriches’ the ISM, that is, it increases the abundance of heavy elements in the ISM. Inflows of gas from the intergalactic medium (IGM) and outflows driven into the IGM by, e.g., supernova explosions, further affect the abundance of heavy elements in the ISM at any given time, typically lowering it. But overall, the evolution is such that the enrichment of the ISM by stellar nucleosynthesis leads to a steady increase in the abundance of heavy elements in the ISM over time. This cycle is shown in the following illustration:
(Credit: Richard Longland)
The previous paragraph presents an overly simplified picture, but is nevertheless a useful starting point to consider the chemical evolution of galaxies. In detail, the dark matter distribution may evolve during the evolution of galaxies and it may itself be affected by the outflows driven by supernovae and other feedback processes (see Chapter 7.7). Star-forming gas is observed to form a thin layer in present-day galaxies, but may not have been as thin earlier in the evolution of the Universe. Infalling gas from the IGM typically has a lower heavy-element abundance than the ISM, but galactic fountains (where gas driven out of one part of the galaxy by outflows returns to a different part of the galaxy) can lead to inflows with higher abundances. For some changes in the rate of star formation and the rate at which gas is lost to outflows, the heavy-element abundance of the ISM can even decline in time. And as the previous chapters have made clear, stars are not static, but orbit within the galaxy over their lifetimes. Thus, enriched ejecta may be returned at very different locations from where a star was born.
Observational galactic astronomy is limited to observing the properties of celestial objects at a single moment. For the Milky Way, this is the present (the light travel time from the furthest reaches of the Galactic halo is \(\lesssim 100,000\) years, and thus smaller than the smallest timescale over which galactic evolution operates). We can therefore only directly determine the heavy-metal abundance of the ISM in the Milky Way as it is today (and even this is difficult). External galaxies observed at moderate or high redshift can give a direct view of the ISM abundance at earlier times, but again only at a single time for a given galaxy. We cannot for a single galaxy directly determine how the ISM abundance changes over time. However, while nuclear burning changes the elemental abundances in the inner parts of a star, the abundances of elements in the stellar photosphere near the surface remain by and large unchanged over the lifetime of a star. Because the photospheric abundances are those that are most easily determined, using spectroscopic observations, the observed heavy-metal abundances of stars today therefore measures the abundances of the ISM when the star was formed. Thus, the elemental abundances of a sample of stars formed over the lifetime of a galaxy allows us to indirectly constrain the chemical evolution of the ISM. Such elemental abundances of a sample of stars can be obtained by measuring abundances for individual stars in the Milky Way through, typically, high-resolution spectroscopy or using stellar-population-synthesis modeling of the integrated spectrum of an external galaxy.
As an example of a measurement of the abundancs of different elements, we can look at the abundances of the Sun as determined by high-resolution spectroscopy of the solar photosphere by Asplund et al. (2009), of which only the part up to atomic number \(Z = 50\) is shown in the Figure below:
(Credit: Bergemann & Serenelli 2014)
The abundances by number of atoms of different elements are shown here on a logarithmic scale where the number of hydrogen atoms in the Sun is \(10^{12}\). To obtain the abundances by mass one would have to multiply by the atomic weight of each element (see Equation \ref{eq-chemev-abundance-weight-from-abundance-number}). It is clear that the pattern of abundances is complex, with an overall decline towards higher \(Z\), an alternating even/odd-\(Z\) pattern, a severe dip for the light elements lithium (Li), beryllium (Be), and boron (B) and a peak around iron (Fe). These trends reflect the relative ease with which these elements are created and destroyed in the Universe (Burbidge et al. 1957). The lightest elements in the Universe, He and Li, were created in the early Universe in a process called Big Bang nucleosynthesis. All elements heavier than boron that naturally occur in the Universe are believed to form in stars or stellar explosions, including some of the most extreme and energetic processes in the evolution of stars (beryllium and boron form primarily by the breaking up of heavier elements in the ISM when they collide with cosmic rays, a process called cosmic ray spallation). The processes that create carbon, nitrogen, oxygen, and heavier elements fall in four broad classes:
Type II supernovae (type II SNe): massive stars (\(M \gtrsim 8\,M_\odot\)) explode as core-collapse supernovae when nuclear burning in the star has reached it endpoint through the formation of an iron/nickel core that cannot uphold the gravity of the star. The ensuing explosion returns material enriched through quiescent stellar nucleosynthesis before the explosion (producing most of the abundance of even-\(Z\) elements lighter than silicon) as well material enriched through explosive nucleosynthesis during the explosion (producing odd-\(Z\) elements and even-\(Z\) elements starting with silicon). Type II supernovae enrichment is responsible for much of the abundance of elements between carbon and scandium and some of the abundance of elements from titanium to zinc (Andrews et al. 2017). Overall, about 90% of all of metals in the Universe are produced by type II supernovae.
Type Ia supernovae (type Ia SNe): low-mass stars end their lives as white dwarfs, which are normally stable. But they can become unstable if they accrete mass from a giant companion (now considered to not significantly contribute to the number of type Ia SNe) or when they slowly (through binary inspiraling) or quickly (through direct collisions) merge with another white dwarf, at which point they explode. Heavy metals returned to the ISM by type Ia supernovae are almost all formed in the explosive nucleosynthesis stage, producing elements from silicon to zinc like for type II SNe. The combination of type II and type Ia SNe enrichment accounts for about 99% of all of the metals in the present-day Universe and for all of the abundances of elements from oxygen to zinc.
The s process: Elements heavier than zinc are not produced during the explosive nucleosynthesis occurring during a supernova explosion. Instead, these elements are believed to be produced through neutron-capture processes, where repeated cycles of neutron capture(s) and beta decay allow elements with higher and higher numbers of nucleons to be formed. The way in which this process occurs is classified into two classes depending on whether the neutron captures occur slowly (for the s [slow] process) or rapidly (for the r [rapid] process) compared to the typical timescale for beta decay of the element in question (which can be minutes to thousands of years or more depending on the isotope in question). Both of these broad processes have multiple basic mechanisms contributing them, with the s process occurring during helium-shell burning in asymptotic giant branch (AGB) stars as well as during helium burning in the cores of very massive stars; both of these occuring on timescales of hundreds of thousands of years. The products from s-process enrichment are ejected into the ISM through stellar winds or during the final demise of the star. Besides creating about half of the abundance of elements heavier than zinc, s-process enrichment is also a major contributor to the abundance of carbon and nitrogen in the Universe. However, at \(Z \lesssim Z_\odot/10\), certain AGB stars actually destroy carbon and oxygen and this leads to a net decrease in the carbon and oxygen abundance (Karakas 2010).
The r process: Like the s process, the r process creates elements heavier than iron, but in environments where neutron captures occur on much shorter timescales than that of the relevant beta decay. Such environments must be neutron rich and the sites of the r process still remain somewhat mysterious. Neutron-star mergers were long suspected of playing a role in the r process and this role was confirmed during the first observed neutron-star merger, known as GW170817, whose lightcurve indicated the presence of ejecta rich in elements produced by the r process (Drout et al. 2017; Kasen et al. 2017). However, it remains unclear what the other sites of the r process are, because the observed chemical evolution of r-process elements indicates that a mechanism that returns r-process enriched ejecta quickly after star formation must be present, while neutron-star mergers typically occur only billions of years after star formation. Like the s process, r-process enrichment creates about half of the abundance of heavier than zinc.
The relative contribution to the abundance of elements from these four classes (six when including Big Bang nucleosynthesis and cosmic ray spallation) for a star like the Sun is summarized in the chart below:
(Credit: Jennifer Johnson; CC-BY-SA)
The four classes discussed above each enrich the ISM on different typical times after star formation creates a new population of stars. Some, especially the r process, may consist of sub-classes that have themselves different typical times. Because enrichment due to type II supernovae happens at the end of massive stars’ lives, the timescale is that of the lifetime of high-mass stars, which is about 5 to 40 Myr for the relevant stars (\(8 \lesssim M/M_\odot \lesssim 50\)); this is short compared to other timescales relevant to galactic evolution. Type Ia supernovae occur when a white dwarf formed at the end of the life of a low-mass star explodes and therefore requires at least the lifetime of the most massive stars that form white dwarfs (tens of millions of years) to start occurring. All proposed mechanisms for how stable white dwarfs turn unstable have a typical timescale of a billions year, so type Ia enrichment occurs on Gyr timescales (but with a significant prompt [\(< 0.5\) Gyr] contribution and a wide distribution of possible times; Maoz et al. 2012). The main contribution to the s process comes from AGB stars, a late stage in the stellar evolution of \(M \lesssim 8\,M_\odot\) stars that typically lasts about 1 Myr, but only occurs near the end of a star’s life and AGB s-process enrichment therefore occurs on the wide range of timescales corresponding to the lifetimes of low- and intermediate-mass stars. As discussed already above, r-process enrichment from neutron-star mergers has a typical timescale of Gyr because of the long time required for neutron stars to spiral in towards a merger. But r-process enrichment due to, e.g., collapsing, rotating massive stars would occur on the same few Myr timescale as type II supernovae (Siegel et al. 2019). Because of these different timescales, the picture presented above on the relative contribution from the different classes to the abundance of different elements changes when we look at stars formed at different times, like early after the formation of a galaxy.
The amount of a particular element returned to the ISM by the various enrichment processes is called the yield, defined as the amount of mass of element X per solar mass of stars formed. These yields depend on the mass of the star, its own elemental abundance mix, and possibly on properties such as the initial rotation, mass transfers from binary companions, etc.; yields are calculated using stellar models including all known nucleosynthetic pathways. For example, a core-collapse supernova of a \(35\,M_\odot\) star produces \(\approx5\,M_\odot\) of oxygen and \(\approx1\,M_\odot\) of carbon (Chieffi & Limongi 2004). For studying galactic chemical evolution, what we care about is not the yield of each individual star, but the yield from an entire population of stars formed at the same time (a star-formation event) at a given time after time star formation. To obtain these population yields, we need to integrate the individual-star yields over the distribution of stellar masses (the initial mass function of stars, IMF). As an example, integrating over an entire stellar population, we have that a core-collapse supernova creates \(\approx0.016\,M_\odot\) of oxygen per solar mass of stars formed. This number is much smaller than the yield of an individual high-mass star, because it takes into account the large amount of mass in low-mass stars that do not explode as core-collapse supernovae and therefore return zero oxygen via this process. Population-level yields in detail depend on the ISM abundance of the stellar population at formation. Because the IMF necessary to integrate the individual-star yields into population yields could plausibly change with time and/or environment, yields could also depend on the time at which stars form. Assuming that the IMF is constant in time, the total yield of all metals is \(\approx 0.035\,M_\odot\) for every solar mass of stars formed at solar metallicity, with \(\approx 90\%\) of this contributed by type II supernovae, \(\approx 8\%\) by type Ia supernovae, and \(\approx 1.5\%\) from AGB winds. At \(Z \approx Z_\odot/10\), the total yield is only \(\approx 0.02\,M_\odot\) and this declines further until it bottoms out at \(\approx 0.01\,M_\odot\) at \(Z \approx Z_\odot/100\). The main driver of this decrease towards lower metallicity is the net destruction of carbon and oxygen at low metallicities by AGB stars, because the yields from supernovae are not very sensitive to metallicity except near \(Z\approx Z_\odot\).
The aim of this chapter is to obtain a basic understanding of galactic chemical evolution through simple models that capture much, but not all, of the relevant physics and to use this to comprehend the present-day distribution of metals and the observed ratios of abundances of different elements in a disk galaxy like the Milky Way.
12.2. Classical models of galactic chemical evolution¶
Simple models of galactic chemical evolution consider at any given time for a given region in a galaxy (i) how much gas is accreted from the environment, (ii) what fraction of the gas in the region turns into stars (the star formation efficiency), (iii) how many metals are returned to the interstellar medium from previous generations of stars, and (iv) how much gas is lost from the system. Alternatively, we can prescribe the number of stars formed at any given time—the star formation rate—which in the above is a combination of points (i) and (ii), instead of the gas accretion. In both cases, these four basic processes can be modeled using simple prescriptions. In the classical models that we consider in this section, the gas accretion history is specified.
The classical models simplify the chemical enrichment treatment by considering that enrichment by the different enrichment processes discussed above happens directly after a stellar population forms. This is called the instantaneous recycling approximation. We also assume that any enriched material is immediately mixed into the ISM, which could technically be considered to be part of the instantaneous recycling approximation as well, but we list it explicitly. We further assume that all yields are absolutely constant, that is, they do not depend on time, stellar abundance, or on any other property of the stellar population or galaxy. Then the distinction between the different enrichment processes vanishes and the evolution of the ISM’s elemental abundance becomes the same for all elements. Thus, this approach cannot be used to investigate how abundance ratios of different elements vary over time as observed, but it is useful to get an overall sense of the build-up of heavy elements in galaxies.
While the approximations in the classical models may appear to be strong, they in fact produce realistic trends for a few elements, notably oxygen and magnesium. As we saw above, these elements are produced primarily in type II supernovae, which, because they originate from massive stars, happen so soon after star formation that they can be considered to be instantaneous to a good approximation. Because their abundance in the ISM receives very small contributions from the other enrichment processes, they are therefore well captured by the approximations of the classical models. However, the classical models were, classically, applied to the iron abundance as well. Because the present-day iron abundance in the ISM receives significant contributions from both type II and type Ia SNe enrichment and the latter occurs typically 1 Gyr after star formation, the instantaneous recycling and single-enrichment-process approximations fail to capture significant aspects of the evolution of the iron abundance.
12.2.1. The closed box model¶
The simplest model for galactic chemical evolution considers a small region (or annulus) in a galaxy and removes any inflow or outflow of gas from consideration (e.g., van den Bergh 1958; Salpeter 1959; Schmidt 1959). Thus, all of the gas that turns into stars is assumed to be present initially, this gas turns into stars over time with some efficiency, metals are returned to the ISM through enrichment and mixed into the ISM, and this enriched gas forms the next generations of stars without losing any gas to the environment. Because of this setup, this is called the closed box model (Talbot & Arnett 1971).
Because it is one of the simplest possible galactic chemical evolution models, the closed box model is often used as a baseline to compare more complex models to. It is therefore instructive to derive the evolution of the heavy-metal content in the ISM in the closed box model. For this we consider the evolution of three quantities: the gas mass \(M_g\), the stellar mass \(M_*\), and the gas mass in metals \(M_Z\); the gas mass includes the gas mass in heavy metals and we denote the metallicity as
Because enrichment happens instantaneously, the following sequence of events happens at every time: (i) some gas turns into stars, (ii) all mass that would ever be returned due to stellar evolution and stellar explosions is returned to the ISM, some enriched in heavy-metal content, (iii) all other stellar mass formed is locked up into stars forever. Because of this, at the end of a given time step, any difference in the gas mass from the previous step is mass that has been turned into stars, or
(Note that here and in everything that follows in this chapter, we ignore the fact that the metals returned to the ISM slightly increase its mass; this increase at most a few percent at the highest ISM abundances observed and it therefore produces a negligible effect on the overall ISM mass and its abundances.)
The evolution of \(M_Z\) is slightly more complicated. At any given time, \(M_Z\) decreases as some of the metals in the ISM are incorporated into stars, through \(\dot{M}_{Z,*-} = Z\,\dot{M}_g=-Z\,\dot{M}_*\), while \(M_Z\) increases as metals are yielded through enrichment, as \(\dot{M}_{Z,*+} = p\,\dot{M}_*\). The parameter \(p\) is the population-level yield of the stellar population, the mass in metals returned to the ISM relative to the mass in stars formed for an entire population of stars; above we discussed that \(p \approx0.035\) for the yield of metals of any kind. Note that winds from massive stars also return (“recycle”) mass to the ISM that has the star’s birth metallicity \(Z\), but we will ignore this for now (we will see below that we can incorporate this by a simple redefinition of the yield \(p\)). The full evolution equation of \(M_Z\) is therefore
Combining Equations \eqref{eq-closedbox-dotmg} and \eqref{eq-closedbox-dotmz}, we have that
or in terms of \(Z\)
using that \(\dot{Z} = \dot{M}_Z/M_g - Z\,\dot{M}_g/M_g\). The solution of Equation \eqref{eq-closedbox-dotz-2} is
Thus, we see that the metallicity \(Z\) increases monotonically as the gas mass decreases due to star formation.
To go further and obtain a specific expression for \(Z(t)\) and \(M_g(t)\) we require the fourth ingredient mentioned at the start of this section, a prescription for the fraction of gas turned into stars at each time step, which we have not done so far. Before doing that, however, we can already obtain a few interesting results. The first is an approximate determination of the population-level yield. This we can obtain using the fact that the current gas fraction in the Milky Way is about 10% (see Section 2.2.2), so within the assumptions of the closed box model that means that \(M_g(\mathrm{now})/M_g(0) \approx 0.1\). The present-day abundance of the ISM is approximately solar (Nieva & Przybilla 2012), and \(Z_\odot \approx 0.014\) (Asplund et al. 2009). Then we can determine the yield \(p\) using Equation \eqref{eq-closedbox-zt} and we find that \(p \approx 0.006\). Comparing this to the value that we quoted above, \(p \approx0.035\), we see that \(p\) is substantially lower. As we will discuss below, this is an indication that the closed box model is incorrect and, in particular, that the yield is low due to outflows of gas that are ignored in the closed box framework.
Secondly, we can derive the metallicity distribution in the closed box model without specifying the star formation efficiency. To do this, we note that because the ISM’s metallicity monotonously increases with time, the number of stars with a given metallicity \(Z'\) is given by the number of stars formed at the time that \(Z(t) = Z'\) and, equivalently, the cumulative number of stars up to a given metallicity \(Z'\) is proportional to the amount of mass in stars formed up to this same time. The latter in the closed box model is simply \(M_g(0) - M_g(t)\) and putting this all together we have that
where we have used Equation \eqref{eq-closedbox-zt} in the last step. The metallicity distribution is then
Observationally, stellar metallicities and their distribution are typically represented as \(\log_{10}\left(Z/Z_\odot\right)\), which we will refer to as \([\mathrm{M/H}]\) to make it clear that this refers to the full metallicity rather than the abundance of a specific element. The metallicity distribution in bins of \([\mathrm{M/H}]\) is
where \([\mathrm{M/H}]=\log_{10}\left(Z/Z_\odot\right)\). Note that when we look at the metallicity distribution after a finite time \(t\), we should cut off this distribution at \(Z(t)\) or its \([\mathrm{M/H}]\) equivalent, but to understand the full range of metallicities that could form in the closed box model, we will simply plot the entire distribution for \(t \rightarrow \infty\) below.
Using the yield \(p = 0.006\) obtained above from the present-day mass and abundance of the ISM, we can plot the closed box metallicity distribution:
[7]:
mhs= numpy.linspace(-3.,1.,201)
Zsolar= 0.014
p= 0.006
Zs= Zsolar*10**mhs
mhdist= Zs*numpy.exp(-Zs/p)
# Normalize such that int d mh mhdist = 1
mhdist/= numpy.sum(mhdist)*(mhs[1]-mhs[0])
figsize(7,5)
plot(mhs,mhdist)
xlabel(r'$[\mathrm{M/H}]$')
ylabel(r'$\mathrm{relative\ fraction}$');
The resulting metallicity distribution is broad with a large number of stars below one-tenth solar metallicity and a peak at \([\mathrm{M/H}] = \log_{10}\left(p/Z_\odot \right) \approx -0.368\). We compare this distribution to observations of the metallicity distribution in the solar neighborhood in the section below. If we change the yield \(p\), for example to the value \(p \approx 0.035\) for the enrichment expected from stellar evolution (which would lead to a present-day ISM abundance of \(Z \approx 5.75 Z_\odot\) or \([\mathrm{M/H}] \approx 0.76\)!) or to half of the value that we derived, the metallicity distribution shifts while largely keeping the same shape:
[147]:
mhs= numpy.linspace(-3.,1.,201)
Zsolar= 0.014
def mhdist_closedbox(mhs,p=0.006):
Zs= Zsolar*10**mhs
mhdist= Zs*numpy.exp(-Zs/p)
# Normalize such that int d mh mhdist = 1
mhdist/= numpy.sum(mhdist)*(mhs[1]-mhs[0])
return mhdist
figsize(7,5)
plot(mhs,mhdist_closedbox(mhs,p=0.006),label=r'$p = 0.006$')
plot(mhs,mhdist_closedbox(mhs,p=0.035),label=r'$p = 0.035$')
plot(mhs,mhdist_closedbox(mhs,p=0.003),label=r'$p = 0.003$')
xlabel(r'$[\mathrm{M/H}]$')
ylabel(r'$\mathrm{relative\ fraction}$')
legend(frameon=False,fontsize=18.);
Physically, the reason that many stars formed have low metallicity in the closed box model is as follows. With all of the gas in place at the start of star formation in the closed box model, early star formation and enrichment produces an amount of heavy metals that is easily diluted in the large gas reservoir necessary to form all of the stars that will ever form. Thus, the metallicity of the ISM at early times grows slowly and many stars therefore form from gas that remains low metallicity for a long while. As the gas supply starts to get exhausted at later times, star formation and subsequent enrichment produces the same yield of heavy elements, but they start to be mixed into a smaller and smaller gas supply, such that the metallicity of the gas rapidly increases and fewer and fewer stars form in a given metallicity range. Thus, \(\mathrm{d}N/\mathrm{d}Z\) steadily decreases with \(Z\) (that there is a maximum in \(\mathrm{d}N/\mathrm{d}[\mathrm{M/H}]\) is solely due to the transformation \(Z \rightarrow [\mathrm{M/H}]\)).
That we have been able to determine the metallicity distribution in the closed box model without specifying how much gas turns into stars at each time is a demonstration of a general rule that the metallicity distribution is largely determined by interplay of the gas supply and the enrichment yield, with only a minor dependence on the star-formation history.
To fully solve chemical evolution in the closed box model, we need to specify how efficiently star formation happens. A simple model is that the star formation efficiency timescale \(\tau_*\) is constant over time, where this quantity expresses how fast the available gas is consumed by star formation and it is defined as
(Note that a larger value of \(\tau_*\) corresponds to a lower efficiency of star formation, because it means that it takes longer to consume the available gas). Combining this with Equation \eqref{eq-closedbox-dotmg}, we have that
or
Thus, the gas supply dwindles exponentially for a constant star formation efficiency in the closed box model. The star formation rate \(\dot{M}_*\) similarly declines, using Equation \eqref{eq-closedbox-taustardef}
The evolution of the metallicity \(Z(t)\) is then given by plugging this into Equation \eqref{eq-closedbox-zt}
Thus, the metallicity \(Z\) increases linearly, reaching the peak \(Z=p\) of the \([\mathrm{M/H}]\) distribution after \(t= \tau_*\) and solar metallicity at \(t = \tau_*\,Z_\odot/p\) or \(t \approx 2.33\tau_*\) for the \(p=0.006\) yield that we derive. If we require that the metallicity of the ISM in the solar neighborhood has solar metallicity today at \(t \approx 10\,\mathrm{Gyr}\), we have that \(\tau_* \approx 4.3\,\mathrm{Gyr}\).
For the constant star-formation efficiency in Equation \eqref{eq-closedbox-taustardef}, we can explicity derive \(\mathrm{d} N / \mathrm{d} Z\) as
in agreement with Equation \eqref{eq-closedbox-dndz}.
Alternatively, we could assume an alternative constant star-formation efficiency \(\tilde{\tau}_*\) defined as
which is the previous star formation efficiency \(\tau_*\) multiplied by the ratio of the gas mass to the initial gas mass. Constant \(\tilde{\tau}_*\) therefore corresponds to an increasing \(\tau_*\) with time, or a decreasing star formation efficiency. Then using Equation \eqref{eq-closedbox-dotmg}, we have that
or
or
Thus, the gas mass declines as \(\approx 1/t\) rather than exponentially, as we expected from the decreasing efficiency of star formation with time in this model. The star formation rate in this case is
and the metallicity increases with time as
To obtain a solar-metallicity ISM after 10 Gyr now, we require that \(\tilde{\tau}_* \approx 1.07\,\mathrm{Gyr}\). While the time dependence of the star-formation rate and the ISM’s metallicity is different than in the example above, we still have that
Thus, we explicitly see that for two different assumptions of how gas turns into stars at each time, the metallicity distribution of the formed stars is the same in the closed box model. The following figure shows the different time evolution in these two models:
[120]:
p= 0.006
tau_star_1= 4.3
tau_star_2= 1.07
ts= numpy.linspace(0.,10.,101)
figsize(13,4)
subplot(1,3,1)
plot(ts,numpy.exp(-ts/tau_star_1))
plot(ts,1./(ts/tau_star_2+1))
xlabel(r'$t$')
ylabel(r'$M_g$')
subplot(1,3,2)
plot(ts,1./tau_star_1*numpy.exp(-ts/tau_star_1),label=r'$\tau_* = M_g/\dot{M}_*$')
plot(ts,1./(tau_star_2*(ts/tau_star_2+1)**2.),label=r'$\tau_* = M_g^2/[M_g(0)\,\dot{M}_*]$')
xlabel(r'$t$')
ylabel(r'$\dot{M}_*$')
legend(frameon=False,fontsize=18.)
subplot(1,3,3)
plot(ts,p*ts/tau_star_1)
plot(ts,p*numpy.log(ts/tau_star_2+1))
xlabel(r'$t$')
ylabel(r'$Z$')
tight_layout();
We see that the metallicity increases faster in the second model, because the gas mass declines faster as well. This creates a steeply dropping star-formation rate, such that the resulting metallicity distribution comes out the same.
12.2.2. The G dwarf problem¶
The closed box model was the first proposed model for galactic chemical evolution and one of the first orders of business was to compare its prediction for the metallicity distribution function to the observed distribution for stars near the Sun. This was done in the early 1960s, when stellar spectroscopy was available for only small samples of stars, which were moreover biased in the way they sampled the metallicity distribution (in that stars were often selected for spectroscopy because they seemed to have interesting metal abundances; this remained an issue until the early 2010s!). However, pioneering work by Wallerstein (1962) demonstrated a tight correlation between the metal content of a G-type dwarf (a star similar to the Sun) and its UV excess. The UV excess of a star in this context is the photometric difference between a UV color (e.g., \(U-B\) in the Johnson \(UBV\) system) and that of a typical star with the same optical color (e.g., \(B-V\); Wallerstein 1962 used stars in the Hyades open cluster to define the typical star colors). This correlation is the result of different levels of metal-line blanketing in the UV (e.g., Schwarzschild et al. 1955), with metal-rich stars having so many metal lines in the UV that their UV flux is significantly suppressed with respect to that of more metal-poor stars. Crucially, this relation meant that an unbiased view of the metallicity distribution could be obtained by applying the UV-excess–metallicity correlation to an unbiased photometric sample of stars. Unfortunately, there are no similar relations between photometric observables and abundance ratios (such as [O/Fe]) and obtaining an unbiased determination of the distribution of abundance ratios therefore remained an observational challenge for a long time.
van den Bergh (1962) and Schmidt (1963) were the first to compare the observed distribution of metallicities obtained using the UV-excess method to the predictions from the closed box model. We can re-create what they did by using a more modern sample, the sample of 1,111 FGK stars obtained from the HARPS exoplanet program (Mayor et al. 2003) and analyzed by Adibekyan et al. (2012). To download this catalog, we use the same function that we defined in Chapter 9.2.2 to download a catalog from the Vizier service:
[4]:
# vizier.py: download catalogs from Vizier
import sys
import os, os.path
import shutil
import tempfile
from ftplib import FTP
import subprocess
_ERASESTR= " "
def vizier(cat,filePath,ReadMePath,
catalogname='catalog.dat',readmename='ReadMe'):
"""
NAME:
vizier
PURPOSE:
download a catalog and its associated ReadMe from Vizier
INPUT:
cat - name of the catalog (e.g., 'III/272' for RAVE, or J/A+A/... for journal-specific catalogs)
filePath - path of the file where you want to store the catalog (note: you need to keep the name of the file the same as the catalogname to be able to read the file with astropy.io.ascii)
ReadMePath - path of the file where you want to store the ReadMe file
catalogname= (catalog.dat) name of the catalog on the Vizier server
readmename= (ReadMe) name of the ReadMe file on the Vizier server
OUTPUT:
(nothing, just downloads)
HISTORY:
2016-09-12 - Written as part of gaia_tools - Bovy (UofT)
2017-09-12 - Copied to galdyncourse (yes, same day!) - Bovy (UofT)
2019-12-20 - Copied directly into notes to avoid additional dependency - Bovy (UofT)
"""
_download_file_vizier(cat,filePath,catalogname=catalogname)
_download_file_vizier(cat,ReadMePath,catalogname=readmename)
catfilename= os.path.basename(catalogname)
with open(ReadMePath,'r') as readmefile:
fullreadme= ''.join(readmefile.readlines())
if catfilename.endswith('.gz') and not catfilename in fullreadme:
# Need to gunzip the catalog
try:
subprocess.check_call(['gunzip',filePath])
except subprocess.CalledProcessError as e:
print("Could not unzip catalog %s" % filePath)
raise
return None
def _download_file_vizier(cat,filePath,catalogname='catalog.dat'):
sys.stdout.write('\r'+"Downloading file %s from Vizier ...\r" \
% (os.path.basename(filePath)))
sys.stdout.flush()
try:
# make all intermediate directories
os.makedirs(os.path.dirname(filePath))
except OSError: pass
# Safe way of downloading
downloading= True
interrupted= False
file, tmp_savefilename= tempfile.mkstemp()
os.close(file) #Easier this way
ntries= 1
while downloading:
try:
ftp= FTP('cdsarc.u-strasbg.fr')
ftp.login()
ftp.cwd(os.path.join('pub','cats',cat))
with open(tmp_savefilename,'wb') as savefile:
ftp.retrbinary('RETR %s' % catalogname,savefile.write)
shutil.move(tmp_savefilename,filePath)
downloading= False
if interrupted:
raise KeyboardInterrupt
except:
raise
if not downloading: #Assume KeyboardInterrupt
raise
elif ntries > _MAX_NTRIES:
raise IOError('File %s does not appear to exist on the server ...' % (os.path.basename(filePath)))
finally:
if os.path.exists(tmp_savefilename):
os.remove(tmp_savefilename)
ntries+= 1
sys.stdout.write('\r'+_ERASESTR+'\r')
sys.stdout.flush()
return None
We then write a function to read the catalog:
[5]:
import os, os.path
from astropy.io import ascii
_CACHE_BASEDIR= os.path.join(os.getenv('HOME'),'.galaxiesbook','cache')
_CACHE_VIZIER_DIR= os.path.join(_CACHE_BASEDIR,'vizier')
_ADIBEKYAN_VIZIER_NAME= 'J/A+A/545/A32'
def read_sn_abu(verbose=False):
"""
NAME:
read_sn_abu
PURPOSE:
read the Adibekyan et al. catalog of abundances of
stars in the solar neighborhood
INPUT:
verbose= (False) if True, be verbose
OUTPUT:
pandas dataframe
HISTORY:
2020-10-14 - Written - Bovy (UofT)
"""
# Generate file path and name
tPath= os.path.join(_CACHE_VIZIER_DIR,'cats',
*_ADIBEKYAN_VIZIER_NAME.split('/'))
filePath= os.path.join(tPath,'table4.dat')
readmePath= os.path.join(tPath,'ReadMe')
# download the file
if not os.path.exists(filePath):
vizier(_ADIBEKYAN_VIZIER_NAME,filePath,readmePath,
catalogname='table4.dat',readmename='ReadMe')
# Read with astropy, w/o the .gz
table= ascii.read(filePath,readme=readmePath,format='cds')
return table.to_pandas()
and a function to plot the abundance distribution in differential \(\mathrm{d}N([\mathrm{X/H}]) / \mathrm{d} [\mathrm{X/H}]\) or cumulative \(N(<[\mathrm{X/H}]\)) form:
[6]:
def plot_sn_xhdist(element,plotrange=[-3,1],bins=31):
sndata= read_sn_abu()
hist(sndata['[{}/H]'.format(element.capitalize())],
range=plotrange,bins=bins,histtype='step',
color='k',lw=1.5,
density=True,label=r'$\mathrm{data}$')
xlabel(r'$[\mathrm{{{}/H}}]$'.format(element.capitalize()))
return None
def plot_sn_xhcumuldist(element,plotrange=[-3,1]):
sndata= read_sn_abu()
sort_indx= numpy.argsort(sndata['[{}/H]'.format(element.capitalize())])
plot(sndata['[{}/H]'.format(element.capitalize())][sort_indx],
numpy.arange(1,len(sndata)+1)/len(sndata),
'k-',lw=1.5,label=r'$\mathrm{data}$')
xlabel(r'$[\mathrm{{{}/H}}]$'.format(element.capitalize()))
ylabel(r'$N(<[\mathrm{{{}/H}}])$'.format(element.capitalize()))
return None
Following in van den Bergh (1962) and Schmidt (1963)’s footsteps, we can then compare the distribution of \([\mathrm{Mg/H}]\) to the prediction of the closed box model that we obtained above, using the yield \(p = 0.006\) that we obtained from matching the present-day abundance of the ISM. Note that while the evolution equations that we have derived consider the abundance of all metals \(Z\) and an element like Mg only makes up a small fraction of the total metal content of a star or the ISM, the predicted distribution \(\mathrm{d}N/\mathrm{d}[\mathrm{Mg/H}]\) is equal to \(\mathrm{d}N/\mathrm{d}[\mathrm{M/H}]\). This is because the fact that the abundance ratio \([\mathrm{Mg/H}]\) is normalized to the equivalent ratio in the Sun causes the factor that scales the overall metallicity to the Mg abundance to drop out (this is the case as long as we assume that all elements move in lockstep, as we do in the instantaneous recycling approximation, but it will cease to be true when we consider different enrichment channels; see Section 12.3).
[149]:
figsize(7,5)
plot_sn_xhdist('Mg')
plot(mhs,mhdist_closedbox(mhs,p=0.006),
label=r'$\mathrm{Closed\ box},\ p = 0.006$')
legend(frameon=False,loc='upper left',fontsize=18.);
It is clear that the closed box model predicts a much too wide distribution of metallicities, with the closed box prediction having a significant number of stars with \([\mathrm{Mg/H}] < -1\), while such stars are essentially absent in the data. This was already known to van den Bergh (1962) and Schmidt (1963), was confirmed many times over the ensuing decades (e.g., Pagel & Patchett 1975), and it remains the case even using the larger, unbiased samples of stellar abundances that are now available (e.g., Hayden et al. 2015; Buder et al. 2019).
This excess of metal-poor G-type dwarf stars predicted in the closed box model is known as the G dwarf problem. Because G-type dwarfs are long-lived, with typical lifetimes of 10 Gyr, it is difficult to explain this problem by assuming that old, metal-poor stars have evolved to stellar remnants and subsequent work furthermore demonstrated that the metallicity distribution of K- and M-type dwarfs has the same lack of metal-poor stars (e.g., Mould 1982; Favata et al. 1997; Schlesinger et al. 2012). Thus, the G dwarf problem is a true problem with the closed box model, not an observational effect.
The G dwarf problem is one of those classical problems in astrophysics that has been solved numerous times, but remains presented as a problem (other examples include the missing-satellites problem and the missing baryons problem). This is both because until we have a definitive theory of galaxy formation, the exact resolution of the problem is difficult to isolate, but also because it is a highly educational problem in learning to think about galactic chemical evolution. Many solutions to the G dwarf problem have been proposed over the years.
The reason that the closed box model overpredicts the number of metal-poor stars in the solar neighborhood is because the abundance of metals in the ISM does not rise quickly enough at early times. Therefore, many long-lived stars are formed while the ISM is still metal-poor and these should be present in the present-day sample of G dwarfs. Any solution of the G dwarf problem therefore has to increase the rate of metal production in the past relative to what it is today significantly above the value of this ratio in the closed box model. Early resolutions of the G dwarf problem focused on changing the initial mass function with time, such that more massive stars were produced relative to low-mass stars in the past than there are today (e.g., Schmidt 1963). A larger number of early massive stars increases the metal-production rate in the past, because these massive stars explode as supernovae and enrich the ISM, while the long-lived, low-mass stars that are more commonly formed at later times in this scenario lock up metals from the ISM without enriching the ISM much. However, models of this type that successfully reproduce the observed metallicity distribution are always about to exhaust their gas supply at the present time and combining this with the relatively constant rate of star formation observed in the solar neighborhood, this means that we would have to be living at a very special time in the life of our Galaxy (Thuan et al. 1975). We also now know that the IMF cannot change as significantly over time as required by these models.
Another easy way to resolve the G dwarf problem is to propose that the initial abundance of the ISM at the start of star formation was not the primordial \(Z=0\) value, but that instead the ISM was pre-enriched to some value \(Z_0\) (Truran & Cameron 1971). In this scenario, the equations describing the time evolution of the closed box model remain the same, but the solution of Equation \eqref{eq-closedbox-dotz-2} is now
and following the same steps as those leading up to Equation \eqref{eq-closedbox-dndz}, the cumulative metallicity distribution is now
and zero for \(Z < Z_0\). Comparing this model with \(Z_0 = Z_\odot/10\) with the cumulative metallicity distribution of the data and with the original closed box model, we find (accounting for \(Z_0\) in the estimate of \(p\) based on the present-day gas fraction and the ISM’s present-day solar abundance, only lowers it to \(p = 0.0055\), so we keep using \(p=0.006\) for simplicitly):
[127]:
mhs= numpy.linspace(-3.,1.,201)
Zsolar= 0.014
def mhcumuldist_closedbox(mhs,p=0.006,Z0=Zsolar/10.):
Zs= Zsolar*10**mhs
mhcumuldist= 1.-numpy.exp((Z0-Zs)/p)
mhcumuldist[Zs<Z0]= 0.
# Normalize such that int d mh mhdist = 1
mhcumuldist/= mhcumuldist[-1]
return mhcumuldist
plot(mhs,mhcumuldist_closedbox(mhs,p=0.006,Z0=0.),
label=r"$\mathrm{Closed\ box},\ p = 0.006,\ Z_0 = Z_\odot/10$")
plot(mhs,mhcumuldist_closedbox(mhs,p=0.006,Z0=Zsolar/10.),
label=r"$\mathrm{Closed\ box},\ p = 0.006,\ Z_0 = Z_\odot/10$")
plot(mhs,mhcumuldist_closedbox(mhs,p=0.01,Z0=Zsolar/3.),
label=r"$\mathrm{Closed\ box},\ p = 0.010,\ Z_0 = Z_\odot/3$")
plot_sn_xhcumuldist('Mg')
ylim(0.,1.52)
legend(frameon=False,loc='upper left',fontsize=18.);
We see that abundances now only start at \([\mathrm{Mg/H}] = -1\), but once star formation commences, the abundance of metals increases as fast as it does in the non-pre-enriched closed box model and, therefore, we still end up with far too many metal-poor stars. The reason that the pre-enriched model follows almost the same evolution as the non-pre-enriched model is that in the context of the closed box model, the initial pre-enriched state is almost the same as the state of the non-pre-enriched model once it reaches \(Z_0\). This is because any stars that have formed at \(Z < Z_0\) in the non-pre-enriched model do not affect the evolution of the gas metallicity (because enrichment is instantaneous), so the only difference is in the gas mass at \(Z = Z_0\), which is \(M_g(0)\) in the pre-enriched model and \(M_g(0)\,e^{-Z_0/p}\) in the non-pre-enriched model. However, the total gas mass does not affect the rate at which metals are produced (see Equation \ref{eq-closedbox-dotz-2}). The only real difference is that we need a lower yield \(p\) to get to the present-day ISM abundance, because there is more time to go from \(Z_0\) to \(Z_\odot\) when the gas is pre-enriched to \(Z_0\).
Pre-enrichment therefore does not entirely solve the problem of the over-production of metal-poor stars in the closed box model, but if we are willing to vary the parameters of the model beyond what we have considered so far, it is possible to find a reasonable match to the data. The figure above includes a line that uses a much higher level of initial enrichment \(Z_0 = Z_\odot/3\) (\([\mathrm{Mg/H}] \approx -0.5\)) and a higher yield \(p = 0.01\). As you can see, this provides a good match to the observed cumulative metallicity distribution. In this model, the present-day abundance of the ISM would be \(\approx 2Z_\odot\) or \([\mathrm{Mg/H}] \approx 0.3\), which is uncomfortably high given observational determinations of the abundances of young stars (which should reflect the ISM’s abundance well; Nieva & Przybilla 2012), but not impossible. Thus, if the gas that forms the stars in the solar neighborhood starts out with \(Z \approx Z_\odot/3\), the pre-enriched closed box model matches the observed data quite well.
Historically, the pre-enrichment model was considered a good model, because it was believed that the stellar halo was formed through star formation during the initial collapse of the gas in our Galaxy (Eggen et al. 1962) before this gas settled in a disk and commenced star formation in the disk. In this model, it is then plausible that the gas that settles in the disk is pre-enriched through the star-formation and enrichment cycle that occurred during the formation of the stellar halo. Nowadays, we have much evidence that the stellar halo largely formed through the accretion of small satellite systems, such as dwarf galaxies and globular clusters (Searle & Zinn 1978), with an additional contribution from stars perturbed out of the disk at early times, but no star formation occurring in the stellar halo itself. Pre-enrichment is therefore not a viable model when considering the entire disk. But a more modern version of the pre-enrichment scenario is that there may be different epochs of star formation in the disk, with in particular an early epoch that enriches the gas to \(Z \approx Z_\odot/3\), but whose stars are rare in stellar samples close to the Sun because these old stars are typically closer to the center or far from the midplane where the Sun lies. Fully exploring such a model requires a much more complete model for disk galaxy formation than we are considering here and we will come back to such more complete models in Part IV.
The pre-enrichment scenario is a solution of the G dwarf problem that does not qualitatively change the basic assumptions of the model. Below, we consider possible solutions that do change the assumptions, in particular stepping away from the ‘closed’ nature of the system.
12.2.3. The leaky box model¶
The basic assumption of the closed box model which gives it its name—that the system is closed—has long been known to be manifestly untrue. Galaxies as a whole and galactic neighborhoods within galaxies typically accrete gas from their environment throughout cosmic history and gas gets removed by winds driven by massive stars, supernova explosions, and active galactic nuclei. To resolve the G dwarf problem, it therefore makes sense to consider the effect of inflows and outflows of gas on the closed box model. We start by in this section looking at the effect of outflows of gas, in what is known as the leaky box model.
That galaxies expel a significant amount of gas to their environment is observed through the ubiquitous existence of superwinds and outflows at low (Heckman et al. 1990) and high redshift (Shapley et al. 2003; Weiner et al. 2009) and from the fact that the gas in the halos of nearby disk galaxies is rich in metals and must therefore consist at least partly of enriched gas driven from the disk into the halo (e.g., Peeples et al. 2014). Simulations of galaxy formation in the expanding Universe also consistently show that an amount of gas similar to that consumed by star formation is expelled from disk galaxies through outflows (e.g., Oppenheimer & Davé 2006; Finlator & Davé 2008).
To open up the closed box model to its leaky cousin, the amount of gas lost to outflows is typically parameterized as a fraction \(\eta\), called the outflow efficiency or the mass loading factor, of the gas turned into stars
This parameterization makes sense for all but the largest galaxies, because outflows are driven by the winds or supernova explosions of massive stars, the latter of which occur in the gas-rich environments of star-forming regions. Type Ia supernovae cause less gas to be driven out of the galaxy, as they typically occur long after a star has separated from its birth cloud and therefore have less gas in their surroundings to drive out. In large galaxies (with dark-matter halos well above the Milky Way’s mass of \(10^{12}\,M_\odot\)), significant outflows are driven by active galactic nuclei and these outflows are not directly related to star formation in the galaxy.
To determine the effect of gas outflows on the predictions of the simple box model that we have been considering, we can assume that the outflow efficiency \(\eta\) is constant. Then Equation \eqref{eq-closedbox-dotmg} becomes
The evolution of the amount of mass in metals \(M_Z\) is also affected, because now some amount of gas with metallicity \(Z\) is lost to outflows and Equation \eqref{eq-closedbox-dotmz} becomes
Combining Equations \eqref{eq-leakybox-dotmg} and \eqref{eq-leakybox-dotmz}, we now have that
or in terms of \(Z\)
where \(p' = p/(1+\eta)\). The solution of Equation \eqref{eq-leakybox-dotz-2} is
Thus, as in the closed box model, the metallicity \(Z\) increases monotonically as the gas mass decreases due to star formation and outflows. In fact, comparing Equation \eqref{eq-leakybox-zt} to the equivalent Equation \eqref{eq-closedbox-zt}, we see that the evolution of \(Z\) in the leaky box is the same as that in the closed box if we replace the yield \(p\) with the effective yield \(p' = p/(1+\eta)\). The effective yield is lower than the pure population-level yield of heavy metals returned to the ISM by stars because of the existence of outflows, which remove enriched gas from the galaxy.
As we discussed in the context of the closed box model, winds from massive stars also return mass to the ISM at the stars’ birth abundance \(Z\). Assuming that this return also happens instantaneously and that the return is a recycling fraction \(r\) of the amount of mass turned into stars \(\dot{M}_\mathrm{recycle} = -r\,\dot{M}_*\) (with a minus sign to indicate that this is mass returned to the ISM), it is clear that this recycling leads to exactly the same evolution as the leaky box with \(\eta = -r\). A closed-box with recycling fraction \(r\) therefore has an effective yield \(p' = p/(1-r)\), while if we want to take into account both recycling and outflows, we obtain \(p' = p/(1+\eta-r)\). The recycling fraction is not negligible, but in fact around 40% of the mass of the mass of a stellar population is returned to the ISM at its birth metallicity (Weinberg et al. 2017).
Because the leaky-box model is mathematically equivalent to the closed box model for \(p \rightarrow p'\), the predicted metallicity distribution is the same as that of the closed box model. The only difference is in the interpretation of the yield parameter. In the closed box model, the yield parameter only takes into account enrichment from stars and from using the fact that the present-day ISM has approximately solar abundances, we determined that \(p \approx 0.006\). This was a surprisingly low value compared to the expected yield of \(p \approx 0.035\) from our theoretical understanding of stellar and explosive nucleosynthesis in high- and low-mass stars. Outflows provide one way of resolving this discrepancy: the determined yield in the context of the leaky-box model is \(p' = 0.006\), while the nucleosynthetic yield is \(p = 0.035\) and they are related through the outflow parameter as \(p/p' = 1+\eta\) (or as \(p/p' = 1+\eta-r\) in the presence of recycling). We can therefore reconcile these two yields by having \(\eta \approx 5\), that is, if at all times we expel from the solar neighborhood a significantly larger amount of gas than the gas consumed by star formation, then the ISM’s present-day solar abundance is natural given the nucleosynthetic yields. As we will see below, this combination of the nucleosynthetic yield and the outflow efficiency are what sets the abundance of the ISM at late times even in more complex models of galactic chemical evolution and the outflow efficiency therefore sets the metallicity of the ISM in galaxies in the local Universe: lower outflow efficiencies in higher mass galaxies lead to an increasing metallicity as a function of galaxy mass, a relation known as the mass–metallicity relation (Tremonti et al. 2004).
12.2.4. The accreting box model¶
Beside expelling gas in outflows, galaxies also accrete pristine gas throughout their lives from their intergalactic environment. Such inflows are difficult to detect, but can be seen, e.g., as high-velocity gas clouds falling into the Milky Way (Oort 1970). Cosmological simulations demonstrate that much of a galaxy’s gas is accreted slowly in the form of cold streams (Katz et al. 1996; Keres et al. 2005; Dekel et al. 2009) and that this is the main driver of star formation (Schaye et al. 2010). Models of galactic evolution have therefore included infalling gas for a long time (e.g., Larson 1972; Tinsley 1974). Because in the closed and leaky box models, the metallicity of the ISM is directly determined by the ratio of the current amount of gas to the initial amount, it is clear that we should expect the evolution of the metal content of the ISM to be significantly different in models where new gas is (semi-)continuously supplied to the system.
As a simple example of adding accretion to our closed box model and turning it into the accreting box model, let’s assume that the accretion of gas is such that the total amount of gas in the box is constant in time. In this case, Equation \eqref{eq-closedbox-dotmz} for the evolution of the mass of metals \(M_Z\) in the ISM is still valid, but we can no longer use Equation \eqref{eq-closedbox-dotmg} to write this in terms of the evolution of the gas mass \(M_g\). Instead, we convert Equation \eqref{eq-closedbox-dotmz} to an equation for the time evolution of \(Z\):
from \(\dot{Z} = \dot{M}_Z/M_g - Z\,\dot{M}_g/M_g\) where the second term on the right-hand side is now zero, because the gas mass is constant (\(\dot{M}_g = 0\)). We cannot solve this equation without an assumption about how stars turn into gas (through an assumed star-formation efficiency), but we can generally understand how the metallicity evolves by replacing time with the total mass \(M = M_g + M_*\) in the system as the independent variable. Because \(\dot{M} = \dot{M}_*\), we have that \(\mathrm{d} Z / \mathrm{d} M = \dot{Z} / \dot{M}_*\), such that
This is a linear inhomogenous differential equation, which we can solved by combining the general solution for the homogeneous version (obtained by setting \(p=0\) in this case and solved by \(Z = C\,\exp\left(-M/M_g\right)\) with \(C\) an integration constant) with a particular solution (\(Z = p\) in this case), such that
To determine the integration constant, we note that initially \(M = M_g\) and \(Z=0\), such that \(C = -p\,e\) and therefore the solution of Equation \eqref{eq-accretingbox-dzdm} is
After star formation commences, the total mass quickly becomes larger than the gas mass (e.g., in the Milky Way \(M / M_g \approx 10\) at the present time) and when \(M \gg M_g\), we see that \(Z \approx p\). Therefore, the abundance of the ISM reaches a constant value \(Z_\mathrm{eq}\) that represents an equilibrium between enrichment from supernovae and infall of pristine gas and the value of \(Z_\mathrm{eq}\) is directly set by the stellar yield. In more general models of accretion, this latter statement continues to be true, with \(Z_\mathrm{eq}\) being independent of the rate of infall (a conclusion first reached by Larson 1972). An easy way to see why the equilibrium point has \(Z = p\) is to consider what happens to a volume of gas consumed by star formation: when a volume with abundance \(Z\) and mass \(m_g\) is turned into stars, two things happen: (a) metals produced by enrichment processes are returned with a total mass of \(m_Z = p\,m_g\) and (b) the gas supply is replenished by adding \(m_g\) of pristine \(Z=0\) gas. The net result of this is that a gas mass with metallicity \(Z\) is replaced by one of the same mass and metallicity \(Z'= (p\,m_g)/m_g = p\). The metallicity of the ISM therefore increases until it reaches a steady state at \(Z=p\). Remember that, as always in this chapter, we are working under the assumption that \(Z \ll 1\), such that the mass in metals can be ignored in the total gas mass.
To determine the yield \(p\) in the accreting box model for the solar neighborhood, we can again calibrate the model to the Milky-Way’s present-day gas fraction, \(M/M_g \approx 10\), and the ISM’s solar abundance, \(Z \approx Z_\odot = 0.014\). From Equation \eqref{eq-accretingbox-zM}, we then derive \(p \approx 0.014\) (\(p \approx Z_\odot\), because \(M/M_g \gg 1\)). As in the closed box model, this value is smaller than the value expected from stellar evolution, \(p \approx 0.035\), although not as small as the \(p\approx 0.006\) that we obtained in the closed box model.
To derive the metallicity distribution in the accreting box model, we can again start by noting that the number of stars \(N(<Z)\) with metallicity less than a given value \(Z'\) is proportional to the total mass in stars \(M_*\) when the ISM has reached metallicity \(Z'\), but now we use Equation \eqref{eq-accretingbox-zM} to express that mass in terms of \(Z\) and constants:
The metallicity distribution is then
or in terms of \([\mathrm{M/H}] = \log_{10}\left(Z/Z_\odot\right)\)
We see that the metallicity distribution is sharply peaked at \(Z \approx p\), although in this model, the singularity at \(Z=p\) is never reached, because that would require \(M = \infty\).
By adding outflows to the accreting box of the same form \(\dot{M}_\mathrm{outflow} = \eta\,\dot{M}_*\) as in the leaky box model, it is straightforward to show that Equation \eqref{eq-accretingbox-zM} becomes
and the metallicity distribution is
where as before \(p' = p/(1+\eta)\) is the effective yield. To match the present-day solar abundance of the ISM with \(M/M_g = 10\) and \(p=0.035\) we then require \(\eta = 1.5\) (for \(M/M_g =10\gg 1\), essentially \(Z = p/[1+\eta]\), so \(\eta \approx p/Z_\odot-1\)). Thus, compared to the \(\eta \approx 5\) value that we required in the leaky box model without accretion, less gas needs to be driven from the solar neighborhood for the ISM’s current solar abundance to make sense. Taking into account recycling as well is trickier than before, because recycling affects the total mass \(M\) in a different way than outflows. It is left as an exercise to show that with both outflows and recycling
For the equilibrium, this expression does simplify to again replacing \(1+\eta\) with \(1+\eta-r\) and, because the solar neighborhood is close to equilibrium in this model, for \(r=0.4\), the outflow efficiency therefore rises to \(\eta \approx 2\) to match the solar neighborhood data. As discussed in Section 12.2.3, the combination \(p/(1+\eta-r)\) sets the abundance of the ISM at late times, which continues to be the case even if we let the gas mass increase or decrease slightly over time.
We can compare these accreting box models (without recycling) to the solar-neighborhood metallicity distribution. Because the model is so sharply peaked, measurement uncertainties in the abundances are important in the comparison, because any non-zero uncertainties will smooth out the intrinsic sharp peak in the metallicity distribution. Therefore, we smooth the predicted accreting-box distributions with a Gaussian kernel of \(0.1\,\mathrm{dex}\) to account for the typical uncertainty in the stellar abundance measurements:
[230]:
from scipy.ndimage import gaussian_filter1d
def mhdist_accretingbox(mhs,p=0.006,eta=0.,e_mh=0.1):
Zs= Zsolar*10**mhs
mhdist= Zs/(p/(1+eta)-Zs)
mhdist[Zs > p/(1+eta)]= 0.
mhdist= gaussian_filter1d(mhdist,e_mh/(mhs[1]-mhs[0]))
# Normalize such that int d mh mhdist = 1
mhdist/= numpy.sum(mhdist)*(mhs[1]-mhs[0])
return mhdist
plot_sn_xhdist('Mg',bins=51)
# Make sure to avoid [M/H]=0 where there is a singularity
mhs= numpy.linspace(-3.,1.,200)
plot(mhs,mhdist_closedbox(mhs,p=0.006),
label=r'$\mathrm{Closed\ box},\ p = 0.006$')
plot(mhs,mhdist_accretingbox(mhs,p=0.014,e_mh=0.1),
label=r'$\mathrm{Accreting\ box},\ p = 0.014$')
plot(mhs,mhdist_accretingbox(mhs,p=0.035,eta=1.5,e_mh=0.1),
label=r'$\!\!\!\!\!\!\mathrm{Accreting\ leaky\ box},$'
+r'\\\phantom{hah}'+r'$p = 0.035,\ \eta=1.5$')
legend(frameon=False,loc='upper left',fontsize=18.);
The prediction from the non-leaky and leaky versions of the accreting box model are of course the same, because they have the same effective yield \(p'\). We see that the tail of metal-poor stars that was present in the closed and leaky box models is now absent and that the support of the metallicity distribution—the range over which it is non-zero—in the data and the model is much more similar. The shape of the accreting-box metallicity distributions is is similar to that of the data, although it is narrower to a degree that cannot be easily explained by the measurement uncertainties. From our considerations of the closed box model, it should be clear that the low-metallicity tail in the accreting box model could be populated by having a higher gas mass in the past than today, for example, by letting accretion be lower today than it was in the past. Because the Universe was more dense in the past than it is now, this is in fact what we would expect to happen. Populating the high-metallicity tail with more stars is more difficult, because the yield sets an upper limit on the metallicity in the accreting box model. As we will see below, these stars were most likely not formed near the solar neighborhood, but they were born in a region of the Galaxy where the effective yield (once we take into account the effect of outflows as well) is higher than it is locally.
Thus, we see that by considering the effects of outflows and accretion, we can solve two of the major issues with the closed box model as applied to the solar neighborhood:
The lack of observed metal-poor stars (the G dwarf problem) can be explained through slow accretion of pristine gas: because the initial gas reservoir is small in the accreting box scenario compared that in the closed box, the metallicity of the ISM increases much more rapidly as the same number of enriched ejecta get mixed into a small gas reservoir and few metal-poor stars are formed;
The low abundance of the present-day ISM (\(Z \approx Z_\odot=0.014\)) compared to that expected in the closed box model (\(Z \approx 2.3\,p\) where \(p \approx 0.035\)), is explained by a combination of inflows and outflows: outflows lower the effective yield, with stronger outflows leading to bigger suppression (\(p \rightarrow p'=p/[1+\eta]\)) and inflows cause the ISM’s abundance to quickly reach an equilibrium value at \(Z \approx p'\) at which point enrichment by stellar evolution is balanced by inflow of pristine gas from the galactic environment.
A relatively continuous inflow of gas that keeps the star formation rate in the disk of a galaxy steady combined with outflows that are such that approximately the same amount of gas is expelled at any given time as is consumed by star formation therefore suffice to understand the overall chemical evolution of disk galaxies.
12.3. The evolution of abundance ratios¶
In the classical models of galactic chemical evolution that we have considered so far in this chapter, where all recycling is instantaneous and there is a single, time-independent yield, all elements evolve with time in the same way. The only difference between elements is their yield, because the mass of, say, iron yielded by a stellar population per solar mass of stars formed is lower than the mass of, for example, oxygen. The abundance \(Z\) in all of the models above always behaves as \(Z(t) = p\,f(t)\), where \(f(t)\) is some function of time that does not depend on \(p\) (see equations \ref{eq-leakybox-zt} and \ref{eq-accretingbox-zM}). Therefore, if we use the formalism above to consider the evolution of different elements rather than of the total metallicity, the ratio of the abundances of different elements is, e.g.
The ratio of the abundances of different elements is therefore constant with time and given by the ratio of the yields of the different elements. In this picture, the relative abundance of elements in the Sun is therefore also equal to this constant ratio and once we normalize abundance ratios to those observed in the Sun (as we usually do when discussing observations), all logarithmically-expressed abundance ratios, e.g., \([\mathrm{Mg/Fe}]\), are zero.
In reality, however, as we discussed in Section 12.1 above, elements are produced by qualitatively different enrichment processes that act on different timescales and are described by different physics. Moreover, the way each physical process enriches the ISM may depend on the metallicity of the star, its environment, etc. This causes the relative abundance of elements to change with time and the abundance ratios of stars formed at different times to differ from zero. Observations of stars in the solar neighborhood (e.g., the data from Adibekyan et al. (2012) that we used above or the data that we will use below) shows that abundance ratios of elements produced primarily by type II supernovae (like O, Mg, Si, etc.) range from \([\mathrm{Mg/Fe}] \approx -0.1\) to \([\mathrm{Mg/Fe}]\approx0.4\), while for some other elements and in other parts of the Milky Way, ratios as small as \(-1\) and as large as \(+1\) are observed. Inasmuch as these deviations from zero are caused by variations in the yield coming from stellar evolution, they can teach us about the different physical processes occurring in stars of different ages and types. If the abundance-ratio deviations are due to the different time scales of different enrichment processes, changes in the star-formation rate and gas flows over time can affect the observed abundance ratios, as can galactic dynamical effects that can move stars significantly away from their birth sites before they enrich the ISM. In the latter case, abundance ratios can help us understand galactic evolution more generally.
For elements heavier than nitrogen up to zinc, the most important enrichment processes are enrichment due to type II (or core-collapse) and type Ia supernovae. The population yields for most of these elements—especially oxygen, magnesium, and iron-peak elements—are relatively constant with metallicity and can be approximated reasonably well as being absolutely constant. The main driver in the evolution of these elements is therefore the different timescales of type II and type Ia supernovae: type II SNe enrich the ISM approximately instantaneously, because they only occur for massive stars with lifetimes \(\lesssim 40\,\mathrm{Myr}\), while type Ia SNe have a much wider range of ages at which they occur up to the age of the Universe. Because these elements are both abundant in stars and relatively easy to determine observationally, we will focus on describing their evolution here. But the tools we develop can be applied to other processes as well, for example to understand the origin of the \(r\) process.
12.3.1. Enrichment by type Ia supernovae¶
For type II supernovae, the progenitor system—massive stars (Hirata et al. 1987; Smartt 2009)—have been known for a long time even if the exact details of the supernova explosion are not entirely understood (Janka 2012). However, for type Ia supernovae, which are believed to result from thermonuclear explosions of carbon-oxygen white dwarfs, the progenitor system remains more unclear. The basic problem is that typical carbon-oxygen white dwarfs have masses of \(\approx 0.5\) to \(1\,M_\odot\), while such white dwarfs are stable up to the Chandrasekhar mass of \(\approx 1.4\,M_\odot\) (Chandrasekhar 1931; Chandrasekhar 1939). Something therefore needs to cause the mass of a white dwarf to increase to \(\approx 1.4\,M_\odot\) to get it to explode. Presently, the main progenitor systems considered are:
single-degenerate systems composed of a white-dwarf and red-giant binary where mass from the giant accretes through Roche lobe overflow onto the white dwarf until it reaches \(\approx 1.4\,M_\odot\) (Whelan & Iben 1973); and
double-degenerate systems consisting of two white dwarfs that merge either through gravitational inspiral (Iben & Tutukov 1984; Webbink 1984) or through direct collisions in triple stellar systems (Kushnir et al. 2013).
Various lines of evidence, such as the absence of evidence for a giant companion in pre-explosion images or in the lightcurve of nearby type Ia supernovae (e.g., Kelly et al. 2014; Margutti et al. 2014; Lundqvist et al. 2015) and the shape of the delay-time distribution (the number of supernovae occurring at different times after star formation; see below), now indicate that the double-degenerate scenario most likely explains the majority of type Ia supernovae (Maoz et al. 2014). Regardless of which progenitor system turns out to be correct, it is clear that the processes involved in all scenarios take time: type Ia supernovae should be expected to occur on a broad range of timescales, extending up to the multiple Gyr that it can take to accrete enough matter from a companion to reach the Chandrasekhar limit (in the single-degenerate case) or to gravitationally spiral in (in the double-degenerate scenario). The implication of this for studying galactic chemical evolution is that the instantaneous recycling approximation is a bad one for elements that receive significant contributions from type Ia supernova enrichment (mainly the iron peak elements). Such elements have a significantly different time evolution than those that are mainly produced by core-collapse supernovae (such as oxygen and magnesium) and abundance ratios involving such elements in the ISM will therefore not be constant, but evolve over time.
Because most of the metals returned to the ISM by type Ia supernovae are created in the explosion of the white dwarf, the predicted yields of different elements do not depend much on the progenitor system. The most important ingredient for the purpose of galactic evolution is therefore the distribution of times after star formation at which a type Ia supernova system explodes and returns heavy elements to the ISM. This distribution is given by the delay time distribution (DTD), which gives the Ia supernova rate as a function of time following a burst of star formation. Luckily for us, the delay time distribution is directly observationally accessible, by determining the supernova rate using observed supernovae in galaxies and clusters of galaxies and comparing it to the galaxy (cluster) population’s star formation history (Maoz & Mannucci 2012). The figure below shows a compilation of different measurements of the delay time distribution:
(Credit: Maoz & Mannucci 2012)
It is clear from this figure that a significant fraction of type Ia supernovae explode at very large times (\(\gtrsim 5\,\mathrm{Gyr}\)) after star formation. However, while the data at delay times below 1 Gyr are sparse, it is also clear that the rate of explosions is highest at the shortest delay times. The solid curve shows a \(t^{-1}\) decay of the supernova rate, which fits the data well.
Theoretically, we expect a broad distribution of delay times in both the single- and double-degenerate scenarios. In the single-degenerate scenario, the delay time distribution is mostly determined by the requirement that the accretion rate onto the white dwarf has to be just so (\(\dot{M} \approx 10^{-7}\,M_\odot\,\mathrm{yr}^{-1}\); Nomoto 1982) for the white dwarf to stably grow in mass. This leads the delay times to range from a few hundred Myr to a few Gyr, with an exponential decrease at shorter and longer delay times (Greggio 2005; Bours et al. 2013).
In the double-degenerate scenario, a double white-dwarf binary has to form through a complex process of binary evolution that culminates in a common-envelope phase and the loss of the envelope. The long delay time behavior (\(\gtrsim 1\,\mathrm{Gyr}\)) is determined by the time it takes to merge the resulting distribution of white-dwarf binary separations. Generically, the merger time dependence on the initial separation of the binary in either the gravitational-inspiral or direct-collision channel for binary white-dwarf mergers is so steep that it results in a delay time distribution \(\propto t^{-1}\) (Maoz et al. 2014). At shorter delay times, this power-law behavior has to break down, because each white dwarf in the binary’s main sequence progenitor must have a mass \(\lesssim 8\,M_\odot\) and lifetime \(\gtrsim 40\,\mathrm{Myr}\) to avoid exploding as a type II supernova, which therefore sets the absolute minimum time delay. Main-sequence stars with masses down to \(2\) to \(3\,M_\odot\) likely create white dwarfs that can be part of a type Ia supernova progenitor system and we should therefore expect the delay time distribution to fall below the \(t^{-1}\) merger-time behavior up to delay times equal to the main-sequence lifetimes of these low-mass progenitors (\(\approx 1\,\mathrm{Gyr}\) for a \(2\,M_\odot\) star or \(\approx 350\,\mathrm{Myr}\) for a \(3\,M_\odot\) star; Reid et al. 2002). The birth rate of white dwarfs is essentially set by the initial mass function (\(\mathrm{d}N/\mathrm{d}m \propto m^{-2.35}\) at the relevant masses \(m \gtrsim 2\,M_\odot\)) and the main-sequence lifetime (\(t \propto m^{-2.5}\)); for the power-law forms given, this gives a birth rate \(\mathrm{d}N/\mathrm{d}t \propto t^{-0.46} \approx t^{-1/2}\) (Pritchet et al. 2008). Because this falls below the \(t^{-1}\) merger-time behavior, it limits the rate of type Ia supernovae, so the delay time distribution for the double-degenerate scenario roughly looks like: (a) zero up \(40\,\mathrm{Myr}\), (b) \(\propto t^{-1/2}\) from the white-dwarf birth rate up to a time \(t_c\approx 350\,\mathrm{Myr}\) to \(1\,\mathrm{Gyr}\) depending on the lowest-mass main-sequence progenitor that contributes to an explosion, and (c) \(\propto t^{-1}\) from the steep dependence of the merger time on initial separation in the double white-dwarf binary (Maoz et al. 2014).
To be able to solve the non-instantaneous-recycling models for chemical evolution below, we require a SN Ia DTD that is an exponential starting at some minimum time or a sum of such exponentials (Weinberg et al. 2017). Thus, we need to approximate the observed DTD with a sum of a hopefully small number of exponentials. First, we need a good representation of the observed DTD. To determine that, we can use the observed DTD from Maoz et al. (2012), a recent determination based on type Ia supernovae observed with the SDSS-II survey (Frieman et al. 2008) in a parent sample of galaxies for which spectroscopic star-formation histories are available (Tojeiro et al. 2009). For the observed DTD, the total number of supernova explosions after 14 Gyr (the age of the Universe) per unit of stellar mass formed is \(N_\mathrm{SN}/M = 1.3\times10^{-3}\,M_\odot^{-1}\), a number which is uncertain at the \(50\%\) level (Maoz & Mannucci 2012). Let’s compare the data to a model for the DTD \(\Psi(\Delta t) \propto t^{-1.1}\), which close to the best-fitting single power-law, and a model inspired by the expected double-degenerate DTD from the considerations above. The latter model has \(\Psi(t) \propto t^{-1/2}\) between the minimum delay time and the lifetime of the least massive star that contributes to type Ia SNe (we’ll use 570 Myr, which is the approximate lifetime of a \(2.5\,M_\odot\) star) and \(\Psi(t) \propto t^{-1.4}\), a slightly steeper dependence than the naively expected \(t^{-1}\) long-term behavior which fits the data better. For all models, we’ll assume that the minimum delay time is 50 Myr (assuming it takes at least 10 Myr to go from a \(8\,M_\odot\) binary white dwarf pair to an explosion). First, we implement the different DTD models:
[7]:
def dtd_singlepower(dt,p=1.1,NSN_M=1.3e-3,mindt=0.05,maxdt=14.):
"""A single power-law t^-p DTD with total
between 0 and 14 Gyr = NSN_M SN Ia / Mstar"""
norm= (mindt**-(p-1)-maxdt**-(p-1))/(p-1)
out= NSN_M/norm*dt**-p
out[dt<mindt]= 0.
return out
def dtd_brokenpower(dt,p1=0.5,p2=1.4,breakdt=0.57,
NSN_M=1.3e-3,mindt=0.05,maxdt=14.):
"""Broken power-law DTD with dt^-p1 at dt<breakdt
and dt^-p2 at dt>breakdt"""
norm= (mindt**-(p1-1)-breakdt**-(p1-1))/breakdt**-p1/(p1-1)\
+(breakdt**-(p2-1)-maxdt**-(p2-1))/breakdt**-p2/(p2-1)
out= numpy.zeros(len(dt))
out[dt<breakdt]= (dt[dt<breakdt]/breakdt)**-p1
out[dt>=breakdt]= (dt[dt>=breakdt]/breakdt)**-p2
out[dt<mindt]= 0.
return out*NSN_M/norm
def dblexp_approx(dt,hdt1=0.4,hdt2=3.3,
NSN_M=1.3e-3,mindt=0.05,maxdt=14.):
"""A DTD that is a sum of two equally-weighted exponentials"""
out= numpy.exp(-dt/hdt1)/hdt1/(numpy.exp(-mindt/hdt1)-numpy.exp(-maxdt/hdt1))/2.
out+= numpy.exp(-dt/hdt2)/hdt2/(numpy.exp(-mindt/hdt2)-numpy.exp(-maxdt/hdt2))/2.
out[dt<mindt]= 0.
return out*NSN_M
and then we plot them versus the data:
[72]:
from matplotlib.ticker import FuncFormatter
data_dt= numpy.array([0.21,1.41,8.2])
data_dtd= numpy.array([140,25.1,1.83])*1e-5
data_dtd_err= numpy.array([30,6.3,0.42])*1e-5
data_cumdtd_x= numpy.array([0.42,2.4,14.])
data_cumdtd_y= numpy.cumsum(data_dtd*data_cumdtd_x)
data_cumdtd_y_err= numpy.sqrt(numpy.cumsum(data_dtd_err**2.*data_cumdtd_x))
figsize(7,5)
dts= numpy.linspace(0.051,14.,201)
loglog(dts,dtd_singlepower(dts),
label=r'$\propto t^{-1.1}$')
loglog(dts,dtd_brokenpower(dts),
label=r'$\propto \begin{cases} t^{-0.5}\,,& t < 0.57\,\mathrm{Gyr}\\ t^{-1.4}\,,& t \geq 0.57\,\mathrm{Gyr}\end{cases}$')
loglog(dts,dblexp_approx(dts),
label=r'$\mathrm{exp.\ approx.}$')
errorbar(data_dt,data_dtd,yerr=data_dtd_err,marker='o',color='k',
ls='none')
xlabel(r'$\Delta t\,(\mathrm{Gyr})$')
ylabel(r'$\Psi(\Delta t)\,(\mathrm{SNe\ Ia}\ \mathrm{Gyr}^{-1}\,M_\odot^{-1})$')
gca().xaxis.set_major_formatter(FuncFormatter(
lambda y,pos: (r'${{:.{:1d}f}}$'.format(int(\
numpy.maximum(-numpy.log10(y),0)))).format(y)))
legend(frameon=False,fontsize=18.,loc='lower left');
We see that the broken power-law model fits better, because the data indicate a steeper drop at \(t > 1\,\mathrm{Gyr}\) than would be expected from a single power law. For the effect of type Ia supernova enrichment on galactic chemical evolution, the cumulative number of explosions is important, so we compare the different models to the data as cumulative functions as well:
[93]:
data_cumdtd_x= numpy.array([0.42,2.4,14.])
# integration interval for each observation
intdt= numpy.hstack((data_cumdtd_x[0],numpy.diff(data_cumdtd_x)))
data_cumdtd_y= numpy.cumsum(data_dtd*intdt)
data_cumdtd_y_err= numpy.sqrt(numpy.cumsum(data_dtd_err**2.*intdt**2.))
figsize(7,5)
dts= numpy.linspace(0.0,14.,1001)
d_dt= dts[1]-dts[0]
plot(dts,numpy.cumsum(dtd_singlepower(dts))*d_dt,
label=r'$\propto t^{-1.1}$')
plot(dts,numpy.cumsum(dtd_brokenpower(dts))*d_dt,
label=r'$\propto \begin{cases} t^{-0.5}\,,& t < 0.57\,\mathrm{Gyr}\\ t^{-1.4}\,,& t \geq 0.57\,\mathrm{Gyr}\end{cases}$')
plot(dts,numpy.cumsum(dblexp_approx(dts))*d_dt,
label=r'$\mathrm{exp.\ approx.}$')
errorbar(data_cumdtd_x,data_cumdtd_y,yerr=data_cumdtd_y_err,marker='o',color='k',
ls='none')
axhline(1.3e-3/2.,color='0.7',ls='--')
ylim(-0.000095,0.001495)
xlabel(r'$\Delta t\,(\mathrm{Gyr})$')
ylabel(r'$\int\,\mathrm{d}\Delta t\Psi(\Delta t)\,(\mathrm{SNe\ Ia}\ M_\odot^{-1})$')
legend(frameon=False,fontsize=18.,loc='lower right');
This again shows that there are more late-time explosions in the single power-law model than in the physically-motivated model. The single-power model is in slight tension with the data at \(\Delta t \approx 2\,\mathrm{Gyr}\). The dashed line indicates half of the supernova explosions that will happen; all models agree that these happen within 1 Gyr.
Because the broken power-law model compares better to the data, we will use that in the following sections. The figures above also show the best approximation of the broken power-law DTD as a sum of equally-weighted exponentials (there is no need for the exponentials to have the same contribution to the total number of supernova explosions, but it is simply convenient for the expressions below). This approximation has exponentials with decay timescales \(\tau = 0.4\,\mathrm{Gyr}\) and \(3.3\,\mathrm{Gyr}\).
Armed with the SN Ia DTD, we can model how SN Ia yields are distributed in time. For example, a single type Ia supernova produces \(m_\mathrm{Fe}^{\mathrm{Ia}} \approx 0.77\,M_\odot\) of iron (Iwamoto et al. 1999). To get the differential or cumulative yield of iron from type Ia supernovae, we can therefore simply multiply the curves above by \(m_\mathrm{Fe}^{\mathrm{Ia}}\). The total yield of iron 14 Gyr after star formation just from SNe Ia is then \(p_\mathrm{Fe}^{\mathrm{Ia}} \approx 0.001\). This is about the same as the total iron yield from core-collapse supernovae, which is \(p_\mathrm{Fe}^{\mathrm{II}} \approx 0.0011\) (Chieffi & Limongi 2004). The long term yield of iron from different supernova channels is therefore the same, but all of the iron enrichment by core-collapse supernovae happens within 40 Myr, while the type Ia supernova enrichment occurs at 40 Myr and beyond, with about half of this part coming at \(\Delta t \gtrsim 1\,\mathrm{Gyr}\).
12.3.2. The delayed closed box model¶
To get started with understanding the impact of the delay in enrichment resulting from type Ia supernova explosions, we can re-visit the closed box model from Section 12.2.1 taking enrichment delays into account. We will focus on understanding the evolution of the iron abundance in the ISM and the ratio of the abundance of a pure core-collapse element, such as oxygen or magnesium, to the iron abundance. For the core-collapse contribution, we will continue to use the instantaneous recycling approximation, that is, enrichment from core-collapse supernovae is assumed to happen immediately after star formation.
The evolution of the pure core-collapse element remains that of the classical closed box model, expressed in terms of the gas mass in Equation \eqref{eq-closedbox-zt}. To obtain the equation that describes the evolution of the iron mass in the ISM, we need to edit Equation \eqref{eq-closedbox-dotz}. Iron enrichment occurs both through core-collapse supernovae, with a yield \(p_\mathrm{Fe}^\mathrm{II}\), and through type Ia supernovae explosions, with a yield that depends on the past star formation rate \(\dot{M}_*\). We will split the contributions of iron to the ISM into that coming from core-collapse supernovae and that from type Ia supernovae and we will track them separately. Therefore, we will consider the evolution of \(M_\mathrm{Fe}^{\mathrm{II}}\) and \(M_\mathrm{Fe}^{\mathrm{Ia}}\), the iron mass contributed by core-collapse supernovae and by type Ia, respectively. And similar we will consider the abundances \(Z_\mathrm{Fe}^{\mathrm{II}}\) and \(Z_\mathrm{Fe}^{\mathrm{Ia}}\). While in reality we can of course not know which type of supernova produced a given iron atom, for the purpose of the calculation it is useful to imagine that we can (e.g., pretend that different types of supernovae give rise to different isotopes of iron) and that we add them together in the end to obtain the total iron abundance. The evolution of \(Z_\mathrm{Fe}^{\mathrm{II}}\) is the same as that of the pure core-collapse element considered in the classical closed box model, using the yield \(p_\mathrm{Fe}^\mathrm{II}\).
To solve for the evolution of \(Z_\mathrm{Fe}^{\mathrm{Ia}}\), we need to include iron enrichment from all past episodes of star formation at a given time and add it to the evolution equation \eqref{eq-closedbox-dotz}. At a time \(\Delta t\) after an instantaneous burst of star formation, \(m_\mathrm{Fe}^{\mathrm{Ia}}\,\Psi(\Delta t)\) per unit time of iron is added to the ISM from type Ia supernovae explosions (where \(m_\mathrm{Fe}^{\mathrm{Ia}}\) is the mass of iron produced by a single SN Ia and \(\Psi(\Delta t)\) is the delay time distribution). To obtain the total yield of iron from past star formation at time \(t\), we therefore add up the yields from all past star formation. Using the notation from Section 12.2.1, we write this as
We can write this in terms of the total yield \(p_\mathrm{Fe}^{\mathrm{Ia}}\) from type Ia supernovae. We defined this above as \(p_\mathrm{Fe}^{\mathrm{Ia}} = m_\mathrm{Fe}^{\mathrm{Ia}}\,\int_0^{14\,\mathrm{Gyr}}\mathrm{d}\Delta t\,\Psi(\Delta t)\), but because we will use exponential DTDs that fall off steeply beyond 14 Gyr, we can approximate this by extending the integral to infinity, in which case Equation \eqref{eq-chemev-delayclosed-firstdMZFe} can be written as
where the DTD-averaged star formation rate \(\langle \dot{M}_*\rangle(t)\) is defined as
This DTD-averaged star-formation history replaces the instantaneous star-formation history in Equation \eqref{eq-closedbox-dotz}. If the DTD were a delta function at \(t=0\), that is, for instantaneous recycling, \(\langle \dot{M}_*\rangle(t) = \dot{M}_*(t)\). If the DTD were a shifted delta function at some delay \(T\), then we would find \(\langle \dot{M}_*\rangle(t) = \dot{M}_*(t-T)\). That is, the iron enrichment at time \(t\) would be the yield times the amount of mass turned into stars at a time \(T\) in the past rather than today. For realistic type Ia DTDs, the DTD-averaged star-formation history gives the effective mass of stars formed in the past for the purpose of considering iron enrichment of the ISM, such that the iron enrichment can be written in the simple form of Equation \eqref{eq-chemev-delayclosed-seconddMZFe}.
The full evolution equation for \(M_{Z_\mathrm{Fe}^{\mathrm{Ia}}}\) is obtained by adding the decrease in \(M_{Z_\mathrm{Fe}^{\mathrm{Ia}}}\) from the metals incorporated into stars, which similar to before is simply \(\dot{M}_{Z_\mathrm{Fe}^{\mathrm{Ia}},*-} = -Z_\mathrm{Fe}^{\mathrm{Ia}}\,\dot{M}_*\) and therefore
To make progress in solving this equation, we need to know the star-formation history \(\dot{M}_*\) so we can compute the DTD-averaged star formation rate. As we discussed in Section 12.2.1, the star formation rate is not uniquely determined unless we specify how efficiently gas turns into stars. To be able to solve for the evolution of \(Z_\mathrm{Fe}^{\mathrm{Ia}}\) analytically, the model from Section 12.2.1 with a constant star formation efficiency time scale \(\tau_*\) is convenient. We showed previously that this model, which has constant
gas mass evolution that is
and a star-formation history
We further assume that the DTD is exponential
where \(\Delta t^{\mathrm{Ia}}_\mathrm{min}= 50\,\mathrm{Myr}\) is the minimum delay time. The DTD-averaged star-formation history is then
where we have introduced \(\tilde{\tau}\) defined through \(\tilde{\tau}^{-1} = \tau_{\mathrm{Ia}}^{-1}-\tau_*^{-1}\) and the resulting expression holds for \(t > \Delta t^{\mathrm{Ia}}_{\mathrm{min}}\) (with \(\langle \dot{M}_*\rangle(t)=0\) at earlier times). Plugging this into Equation \eqref{eq-chemev-delayclosed-dmZfull-1}, we find that
Following the same steps as in going from Equation \eqref{eq-closedbox-dotz} to Equation \eqref{eq-closedbox-dotz-2}, we then have the evolution equation for \(\dot{Z}_\mathrm{Fe}^{\mathrm{Ia}}\)
We can solve this equation simply through integrating both sides from \(t'=0\) to \(t'=t\), using that \(\dot{M}_g/ M_g = -1/\tau_*\) and that \(Z_\mathrm{Fe}^{\mathrm{Ia}} = 0\) up to \(t=\Delta t^{\mathrm{Ia}}_{\mathrm{min}}\):
at \(t > \Delta t^{\mathrm{Ia}}_{\mathrm{min}}\) and zero at earlier times. That \(Z_\mathrm{Fe}^{\mathrm{Ia}}\) is zero until \(\Delta t^{\mathrm{Ia}}_{\mathrm{min}}\) makes sense, because it takes that long for the first SNe Ia to explode. Taylor expanding the exponential for \(t-\Delta t^{\mathrm{Ia}}_{\mathrm{min}} \ll \tilde{\tau}\) shows that \(Z_\mathrm{Fe}^{\mathrm{Ia}} \propto (t-\Delta t^{\mathrm{Ia}}_{\mathrm{min}})^2\) at small \(t-\Delta t^{\mathrm{Ia}}_{\mathrm{min}}\), so type Ia supernova enrichment turns on slowly.
To apply this model to the solar neighborhood abundance distribution, we need to set the parameters \(\tau_*\) and \(\tau_{\mathrm{Ia}}\). For simplicity, we will for now represent the DTD using a single exponential with \(\tau_{\mathrm{Ia}} = 1.5\,\mathrm{Gyr}\). To set \(\tau_*\), we can use the fact that in the context of the closed box model, the present-day ratio of stellar to gas mass is equal to \(M_g(t)/M_g(0) = e^{-t/\tau_*}\) and the fact that this ratio is about 0.1 therefore requires \(\tau_* \approx 4.5\,\mathrm{Gyr}\) (note that this value only applies in the closed box model). We can then compare the evolution of \(Z_\mathrm{Fe}^{\mathrm{Ia}}\) to that of \(Z_\mathrm{Fe}^{\mathrm{II}}\), which is given by Equation \eqref{eq-chemev-Zt-closedbox-consttauSFE} with \(p = p_\mathrm{Fe}^{\mathrm{II}} \approx 0.0011\):
[15]:
def Z_Fe_Ia_closedbox(t,p_Fe_Ia=0.001,tau_Ia=1.5,tau_star=4.5,mindt=0.05):
tau_tilde= 1./(1./tau_Ia-1./tau_star)
out= p_Fe_Ia*tau_tilde/tau_Ia/tau_star\
*numpy.exp(mindt/tau_star)\
*(t-mindt-tau_tilde*(1.-numpy.exp(-(t-mindt)/tau_tilde)))
out[t<mindt]= 0.
return out
def Z_Fe_II_closedbox(t,p_Fe_II=0.0011,tau_star=4.5):
return p_Fe_II*t/tau_star
ts= numpy.linspace(0.,10.,201)
figsize(7,5)
plot(ts,Z_Fe_Ia_closedbox(ts),label=r'$\mathrm{Ia\ SNe}$')
plot(ts,Z_Fe_II_closedbox(ts),label=r'$\mathrm{II\ SNe}$')
xlabel(r'$t\,(\mathrm{Gyr})$')
ylabel(r'$Z_\mathrm{Fe}$')
legend(frameon=False,fontsize=18.);
In the closed box model, the contribution from instantaneous core-collapse supernovae to the enrichment grows linearly in time, while that from the delayed type Ia supernovae turns on slowly. However, once it turns on, the type Ia enrichment increases faster with time and eventually overtakes the core-collapse contribution. Why and how this happens can be understood by looking at the behavior of Equation \eqref{eq-chemev-delayclosed-dZdt} at \(t \gg \tilde{\tau}\), using the fact that for \(\tau_* \gg \tau_{\mathrm{Ia}}\), \(\tilde{\tau} \approx \tau_{\mathrm{Ia}}\):
Comparing this to Equation \eqref{eq-closedbox-dotz-2} which gives the behavior of \(\dot{Z}_\mathrm{Fe}^{\mathrm{II}}\), we see that the \(t \gg \tilde{\tau}\) behavior is essentially that of the instantaneous recycling approximation, but with the yield given by \(p_\mathrm{Fe}^{\mathrm{Ia}}\,e^{\Delta t^{\mathrm{Ia}}_{\mathrm{min}}/\tau_*}\) (to get quantitative agreement with the slope in the figure above, one needs to account for the \(\tilde{\tau} / \tau_{\mathrm{Ia}}\) factor as well). Physically, this happens because for \(\tau_* \gg \tau_{\mathrm{Ia}}\), all type Ia supernovae explode before the star-formation rate changes significantly and the DTD-averaged star-formation rate is therefore about the same as the instantaneous star-formation rate; thus, this situation is equivalent to an instantaneous recyling model.
In the opposite limit, where \(\tau_* \ll \tau_{\mathrm{Ia}}\) and star formation thus proceeds more rapidly than enrichment from type Ia supernovae, Equation \eqref{eq-chemev-delayclosed-dZdt} at \(t \gg \tilde{\tau}\) becomes
because now \(\tilde{\tau} \approx -\tau_*\). Thus, the effective yield (the part of the expression that multiplies \(\dot{M}_g/M_g\)) increases exponentially with time, leading to exponential growth of the iron abundance \(Z_\mathrm{Fe}^{\mathrm{Ia}} \approx (p_\mathrm{Fe}^{\mathrm{Ia}}\,\tau_*/\tau_{\mathrm{Ia}})\,e^{t/\tau_*}\). This leads to a very high iron abundance; for example, let’s compare the behavior of \(Z^{\mathrm{Ia}}_\mathrm{Fe}\) for \(\tau_*/\tau_{\mathrm{Ia}} = 3\) to that of \(\tau_*/\tau_{\mathrm{Ia}} = 1/3\):
[16]:
ts= numpy.linspace(0.,10.,201)
figsize(7,5)
semilogy(ts,Z_Fe_Ia_closedbox(ts,tau_star=4.5,tau_Ia=1.5),label=r'$\tau_*/\tau_{\mathrm{Ia}} = 3$')
semilogy(ts,Z_Fe_Ia_closedbox(ts,tau_star=0.5,tau_Ia=1.5),label=r'$\tau_*/\tau_{\mathrm{Ia}} = 1/3$')
xlabel(r'$t\,(\mathrm{Gyr})$')
ylabel(r'$Z^{\mathrm{Ia}}_\mathrm{Fe}$')
legend(frameon=False,fontsize=18.);
We have had to switch to logarithmically-scaled axes here because the increase in the iron abundance is so fast for \(\tau_*/\tau_{\mathrm{Ia}} = 1/3\). Of course, this is not a physical model, because the mass in iron of the ISM cannot exceed the total mass of the ISM and therefore \(Z\leq 1\); the high \(Z\) obtained in this model breaks our assumption that \(Z \ll 1\). The physical reason for the exponential increase in the iron abundance when \(\tau_* \ll \tau_{\mathrm{Ia}}\) is that in this case the star formation rate has significantly dropped by the time many of the type Ia supernova explosions happen from earlier star formation. Thus, the DTD-averaged star-formation rate is much higher than the instantaneous star-formation rate, which leads to much higher effective yield than would be expected if enrichment were instantaneous. In the closed box model this effect is amplified by the fact that the gas mass has also decreased significantly between when the stars formed and when they explode as supernovae; the same amount of iron produced by these supernovae then leads to much higher fractional increase of the ISM’s iron abundance.
To confront the predictions from the closed box model with enrichment from both instantaneous core-collapse supernovae and delayed type Ia supernovae, we can calculate the expected relation between the iron abundance \(Z_\mathrm{Fe}\) and the relative abundance of oxygen to that of iron \(Z_\mathrm{O}/Z_\mathrm{Fe}\). For this, we need to compute \(Z_\mathrm{O}\), which assuming that oxygen is only produced by type II supernovae follows Equation \eqref{eq-chemev-Zt-closedbox-consttauSFE} with \(p = p_\mathrm{O}^{\mathrm{II}} \approx 0.016\) (Chieffi & Limongi 2004). To obtain the iron abundance \(Z_\mathrm{Fe}\), we add up contributions from the type Ia and type II enrichment channels: \(Z_\mathrm{Fe} = Z_\mathrm{Fe}^{\mathrm{II}}+Z_\mathrm{Fe}^{\mathrm{Ia}}\). We will also consider the two-exponential approximation to the more realistic type Ia DTD from Section 12.3.1, which splits the DTD into two equal parts with decay timescales \(\tau_{\mathrm{Ia},1}\) and \(\tau_{\mathrm{Ia},2}\). In this case, \(Z_\mathrm{Fe} = Z_\mathrm{Fe}^{\mathrm{II}}+Z_\mathrm{Fe}^{\mathrm{Ia},1}+Z_\mathrm{Fe}^{\mathrm{Ia},2}\), where \(Z_\mathrm{Fe}^{\mathrm{Ia},i}\) is the contribution from type Ia supernovae computed using the \(\tau_{\mathrm{Ia},i}\) exponential. Observationally, abundances and abundance ratios are expressed as logarithmic quantities normalized to the solar abundance, with the additional complication that abundances are weighted by number of atoms rather than by mass. That is, the observed iron abundance \([\mathrm{Fe/H}]\) is defined as
and the oxygen-to-iron ratio \([\mathrm{O/Fe}]\) is
Because the quantity \(Z_\mathrm{Fe}\) that we compute is the mass of iron relative to the total mass rather than the number of iron atoms, the easiest way to convert \(Z_\mathrm{Fe}\) to \([\mathrm{Fe/H}]\) is to express the solar iron abundance as a relative mass of iron in the Sun as well (and similar for \(Z_\mathrm{O}\)). The abundances of different elements in the Sun is typically on the logarithmic scale where the number of hydrogen atoms is \(10^{12}\) (see Section 12.1 above), the abundance by number of element \(X\) is \(x = 12 + \log_{10}\left(\mathrm{X/H}\right)\) and the conversion to a logarithmic mass fraction is
where \(\mu\) is the atomic weight of element \(X\) and \(0.71\) is the mass fraction of hydrogen in the Sun. The atomic weight of iron is 55.85 and the abundance of iron in the Sun is \(x_{\mathrm{Fe}} = 7.5\) while the atomic weight of oxygen is 16 and its abundance in the Sun is \(x_{\mathrm{O}} = 8.69\) (Asplund et al. 2009), so we have that
With these numbers we can then compute \([\mathrm{Fe/H}] = \log_{10} Z_{\mathrm{Fe}} - \log_{10} Z_{\mathrm{Fe},\odot}\) and similar for oxygen. We can then compute the relation between \([\mathrm{O/Fe}]\) and \([\mathrm{Fe/H}]\) (the track) for our two different DTD representations in the delayed closed box model:
[17]:
def feh_closedbox(t,p_Fe_Ia=0.001,p_Fe_II=0.0011,
tau_Ia=0.4,tau_Ia_2=3.3,mindt=0.05,
tau_star=4.5,solar_Fe=7.5):
tau_tilde= 1./(1./tau_Ia-1./tau_star)
out= p_Fe_Ia*tau_tilde/tau_Ia/tau_star\
*numpy.exp(mindt/tau_star)\
*(t-mindt-tau_tilde*(1.-numpy.exp(-(t-mindt)/tau_tilde)))
# Second exponential Ia component
if not tau_Ia_2 is None:
out/= 2.
tau_tilde= 1./(1./tau_Ia_2-1./tau_star)
out+= p_Fe_Ia/2.*tau_tilde/tau_Ia_2/tau_star\
*numpy.exp(mindt/tau_star)\
*(t-mindt-tau_tilde*(1.-numpy.exp(-(t-mindt)/tau_tilde)))
out[t<mindt]= 0.
# Add CCSN contribution
out+= p_Fe_II*t/tau_star
# Convert to [Fe/H]
logZFe_solar= -2.90+solar_Fe-7.50
return numpy.log10(out)-logZFe_solar
def oh_closedbox(t,p_O_II=0.016,
tau_star=4.5,solar_O=8.69):
# Just CCSN contribution
out= p_O_II*t/tau_star
# Convert to [O/H]
logZO_solar= -2.25+solar_O-8.69
return numpy.log10(out)-logZO_solar
ts= numpy.geomspace(1e-3,10.,201)
figsize(7,5)
plot(feh_closedbox(ts),oh_closedbox(ts)-feh_closedbox(ts),
label=r'$\Psi(\Delta t) \propto \begin{cases} t^{-0.5}\,,& t < 0.57\,\mathrm{Gyr}\\ t^{-1.4}\,,& t \geq 0.57\,\mathrm{Gyr}\end{cases}$')
plot(feh_closedbox(ts,tau_Ia=1.5,tau_Ia_2=None),
oh_closedbox(ts)-feh_closedbox(ts,tau_Ia=1.5,tau_Ia_2=None),
zorder=-1,
label=r'$\Psi(\Delta t) \propto e^{-\Delta t / [1.5\,\mathrm{Gyr}]}$')
xlabel(r'$[\mathrm{Fe/H}]$')
ylabel(r'$[\mathrm{O/Fe}]$')
legend(frameon=False,fontsize=18.);
We see that at low \([\mathrm{Fe/H}]\), the relative abundance of oxygen to that of iron is constant and forms a plateau. This occurs, because at early times (low \([\mathrm{Fe/H}]\)), enrichment by type Ia supernovae is negligible and both iron and oxygen are enriched through instantaneous type II supernova enrichment. The relative abundance is then simply equal to the relative ratio of the core-collapse yields for oxygen and iron. However, after about 1 Gyr, enrichment from type Ia supernovae becomes an important contributor to \(Z_{\mathrm{Fe}}\) and the iron abundance therefore starts to increase faster, driving down the relative ratio of oxygen to iron. The exact way in which this happens depends on the DTD, as the figure above shows, although the exact representation of the DTD only makes a small quantitative difference in the location in the track, with the qualitative behavior independent on the details of the DTD. From now on, we will only consider the two exponential approximation to the broken-power-law DTD model.
We can compare the predicted track to data on the abundance ratios of stars in the solar neighborhood. The oxygen abundance is tricky to measure using stellar spectra, because oxygen absorption lines are rare and difficult to interpret in the optical region of the electromagnetic spectrum where most high-resolution spectrographs operate. But we will use a sample of \(\approx 800\) solar-neighborhood stars from Ramirez et al. (2013), which we can load from the Vizier service using the function to download catalogs from Vizier in Section 12.2.2 above. We first write a function to download and process the catalog:
[11]:
import os, os.path
from astropy.io import ascii
from astropy.table import join as table_join
_CACHE_BASEDIR= os.path.join(os.getenv('HOME'),'.galaxiesbook','cache')
_CACHE_VIZIER_DIR= os.path.join(_CACHE_BASEDIR,'vizier')
_RAMIREZ_VIZIER_NAME= 'J/ApJ/764/78'
def read_sn_o_abu(verbose=False):
"""
NAME:
read_sn_o_abu
PURPOSE:
read the Ramirez et al. catalog of oxygen abundances of
stars in the solar neighborhood
INPUT:
verbose= (False) if True, be verbose
OUTPUT:
pandas dataframe
HISTORY:
2020-11-20 - Written - Bovy (UofT)
"""
# Generate file path and name
tPath= os.path.join(_CACHE_VIZIER_DIR,'cats',
*_RAMIREZ_VIZIER_NAME.split('/'))
# First download table with [Fe/H]
filePath_FEH= os.path.join(tPath,'table4.dat')
readmePath= os.path.join(tPath,'ReadMe')
# download the file
if not os.path.exists(filePath_FEH):
vizier(_RAMIREZ_VIZIER_NAME,filePath_FEH,readmePath,
catalogname='table4.dat',readmename='ReadMe')
# Then table with [O/H]
filePath_OH= os.path.join(tPath,'table6.dat')
# download the file
if not os.path.exists(filePath_OH):
vizier(_RAMIREZ_VIZIER_NAME,filePath_OH,readmePath,
catalogname='table6.dat',readmename='ReadMe')
# Read with astropy
table_feh= ascii.read(filePath_FEH,readme=readmePath,format='cds')
table_oh= ascii.read(filePath_OH,readme=readmePath,format='cds')
out= table_join(table_feh,table_oh)
out= out[out['Sample'] == 'main']
return out.to_pandas()
and then a function to plot the \([\mathrm{O/Fe}]\) versus \([\mathrm{Fe/H}]\) distribution of the stars:
[12]:
def plot_sn_ofe_feh(xrange=[-2.05,0.65],yrange=[-0.15,0.6]):
sndata= read_sn_o_abu()
sndata= sndata[(sndata['e_[Fe/H]2'] < 0.15)*(sndata['e_[O/H]N'] < 0.05)]
plot(sndata['[Fe/H]2'],sndata['[O/H]N']-sndata['[Fe/H]2'],
'ko',ms=2.,label=r'$\mathrm{data}$')
xlabel(r'$[\mathrm{Fe/H}]$')
ylabel(r'$[\mathrm{O/Fe}]$')
xlim(*xrange)
ylim(*yrange)
return None
We can then confront the delayed closed box model with \(\tau_* = 4.5\,\mathrm{Gyr}\) with this data:
[18]:
ts= numpy.geomspace(1e-3,10.,201)
figsize(7,5)
plot_sn_ofe_feh()
plot(feh_closedbox(ts),oh_closedbox(ts)-feh_closedbox(ts),lw=2.)
xlabel(r'$[\mathrm{Fe/H}]$')
ylabel(r'$[\mathrm{O/Fe}]$');
We see that while the plateau at low \([\mathrm{Fe/H}]\) is indeed at the upper limit of \([\mathrm{O/Fe}]\) in the data, the drop towards lower \([\mathrm{O/Fe}]\) at higher \([\mathrm{Fe/H}]\) starts at too low \([\mathrm{Fe/H}]\) and is too slow. While the data go down to \([\mathrm{O/Fe}] \approx 0\) and perhaps a bit below (which is expected, because the Sun has to form at some point and the Sun has \([\mathrm{O/Fe}] = 0\)), the delayed closed box model only reaches \([\mathrm{O/Fe}]\approx 0.2\). While the exact values of the plateau and the final \([\mathrm{O/Fe}]\) depend on the yields from the different types of supernovae, which are somewhat uncertain, the discrepancy is beyond that that can be easily explained through theoretical uncertainties and the shallow drop in particular is in strong conflict with the data. Thus, the closed box model can neither explain the observed metallicity distribution (the G dwarf problem, Section 12.2.2) or the track of solar-neighborhood stars in \([\mathrm{O/Fe}]\) versus \([\mathrm{Fe/H}]\)!
12.3.3. The delayed leaky box model¶
We can extend the closed box model with delayed type Ia supernova enrichment with outflows of gas and derive the delayed version of the leaky box model from Section 12.2.3. Again defining the outflow efficiency \(\eta \equiv \dot{M}_{\mathrm{outflow}}/\dot{M}_*\) and assuming it to be constant and combining this now with the constant star-formation-efficiency model (Equation \ref{eq-closedbox-taustardef}; \(\tau_* = M_g/\dot{M}_*\), the evolution of the gas mass becomes
The gas mass therefore evolves as
where \(\tau_{\mathrm{dep}} = \tau_*/(1+\eta)\) is the depletion timescale. The star formation rate follows from \(\tau_* = M_g/\dot{M}_*\)
For this star formation rate, the DTD-averaged star formation rate is still given by the expression in Equation \eqref{eq-chemev-DTDMstar-closedbox}, except that now \(\tilde{\tau}^{-1} = \tau_{\mathrm{Ia}}^{-1}-\tau_{\mathrm{dep}}^{-1}\). The evolution equation for the Ia contribution to iron, \(\dot{M}_{Z_{\mathrm{Fe}}^{\mathrm{Ia}}}\) gains an outflow contribution as in Equation \eqref{eq-leakybox-dotmz-1} and thus becomes
Similar to going from \(\dot{M}_Z\) in Equation \eqref{eq-leakybox-dotz} to \(\dot{Z}\) in Equation \eqref{eq-leakybox-dotz-2} for the classical leaky box model, the evolution equation for \(\dot{Z}_{\mathrm{Fe}}^{\mathrm{Ia}}\) is then the same expression as that in Equation \eqref{eq-chemev-delayclosed-dZdt} but replacing the yield \(p_{\mathrm{Fe}}^{\mathrm{Ia}}\) with the effective yield \(p_{\mathrm{Fe}}^{',\mathrm{Ia}} = p_{\mathrm{Fe}}^{\mathrm{Ia}} / (1+\eta)\) and now \(\dot{M}_g/ M_g = -1/\tau_{\mathrm{dep}}\). The solution of this equation is then also given by the expression in Equation \eqref{eq-chemev-delayclosed-Zt} with the same \(p_\mathrm{Fe}^{\mathrm{Ia}} \rightarrow p_{\mathrm{Fe}}^{',\mathrm{Ia}}\) replacement and additionally replacing \(\tau_* \rightarrow \tau_{\mathrm{dep}} = \tau_*/(1+\eta)\). Thus, similar to the relation between the classical closed and leaky box models, the delayed leaky box model has a similar evolution as the delayed closed box model, but with effective yield and evolution time scale that are lower by \(1/[1+\eta]\) due to outflows. As discussed in Section 12.2.3, we can treat instantaneous recycling of material at the stars’ birth abundance through stellar winds in the same way, leading to an effective yield \(p_{\mathrm{Fe}}^{',\mathrm{Ia}} = p_{\mathrm{Fe}}^{\mathrm{Ia}} / (1+\eta-r)\), where \(r\) is the recycling fraction, which is about 0.4. In this more general case, the depletion timescale is then \(\tau_{\mathrm{dep}} = \tau_*/(1+\eta-r)\).
The same relation between the leaky and closed box models holds for the evolution of the type II supernova contribution to Fe and to that of oxygen, and the evolution of those \(Z_{\mathrm{Fe}}^{\mathrm{II}}\) and \(Z_{\mathrm{O}}^{\mathrm{II}}\) is therefore given by Equation \eqref{eq-chemev-Zt-closedbox-consttauSFE} with yield parameter \(p_{\mathrm{X}}^{\mathrm{II}} \rightarrow p_{\mathrm{X}}^{',\mathrm{II}} = p_{\mathrm{X}}^{\mathrm{II}}/(1+\eta-r)\). To compare the leaky box model with recycling to the data on the relation between \([\mathrm{O/Fe}]\) and \([\mathrm{Fe/H}]\) in the solar neighborhood, we need to set \(\tau_*\) such that \(M_*/M_g \approx 10\) at \(t=10\,\mathrm{Gyr}\); it is left as an exercise to show that this requires \(\tau_{\mathrm{dep}} = t/\ln\left([1+\eta-r]\,M_*/M_g +1\right)\) and therefore \(\tau_* \approx 9\,\mathrm{Gyr}\). This is twice as large as the \(\tau_* = 4.5\,\mathrm{Gyr}\) that we required in the closed box model. That the required \(\tau_*\) to obtain a given ratio of stellar-to-gas mass is larger in the leaky box makes sense, because when some of the gas is lost to outflows, we need to make fewer stars out of the initial gas supply to get to the same ratio and the star formation efficiency can therefore be lower (and its timescale larger). Comparing the leaky box model to the data, we then get
[19]:
ts= numpy.linspace(0.,10.,201)
def feh_leakybox(t,p_Fe_Ia=0.001,p_Fe_II=0.0011,
tau_Ia=0.4,tau_Ia_2=3.3,mindt=0.05,
tau_star=4.5,solar_Fe=7.5,
eta=2.5,r=0.4):
# Equivalent to closed box with p' = p/(1+eta-r)
# and tau_star'= tau_star/(1.+eta-r)
return feh_closedbox(t,p_Fe_Ia=p_Fe_Ia/(1.+eta-r),
p_Fe_II=p_Fe_II/(1.+eta-r),
tau_Ia=tau_Ia,tau_Ia_2=tau_Ia_2,
mindt=mindt,
tau_star=tau_star/(1.+eta-r),
solar_Fe=solar_Fe)
def oh_leakybox(t,p_O_II=0.016,
tau_star=4.5,solar_O=8.69,
eta=2.5,r=0.4):
# Equivalent to closed box with p' = p/(1+eta-r)
# and tau_star'= tau_star/(1.+eta-r)
return oh_closedbox(t,p_O_II=p_O_II/(1.+eta-r),
tau_star=tau_star/(1.+eta-r),
solar_O=solar_O)
ts= numpy.geomspace(1e-3,10.,201)
figsize(7,5)
plot_sn_ofe_feh(xrange=[-2.19,0.65],yrange=[-0.34,0.6])
plot(feh_closedbox(ts),oh_closedbox(ts)-feh_closedbox(ts),
lw=2.,label=r'$\mathrm{closed\ box}; \tau_* = 4.5\,\mathrm{Gyr}$')
plot(feh_leakybox(ts),oh_leakybox(ts)-feh_leakybox(ts),
lw=2.,label=r'$\mathrm{leaky\ box}; \eta=2.5,\ r=0.4,\ \tau_* = 4.5\,\mathrm{Gyr}$')
plot(feh_leakybox(ts,tau_star=9.),
oh_leakybox(ts,tau_star=9.)-feh_leakybox(ts,tau_star=9.),
lw=2.,label=r'$\mathrm{leaky\ box}; \eta=2.5,\ r=0.4,\ \tau_* = 9.0\,\mathrm{Gyr}$')
xlabel(r'$[\mathrm{Fe/H}]$')
ylabel(r'$[\mathrm{O/Fe}]$')
legend(frameon=False,fontsize=18,loc='lower left');
As expected, the leaky box predictions are similar to those of the closed box model, except that they are shifted to lower metallicities because the effective yields are smaller than the true yields by the factor of \(1/(1+\eta-r)\). It is clear that adding outflows to the closed box model does not improve the agreement with the data on the abundance ratios.
12.3.4. The delayed accreting box¶
Next, we can add inflows of gas to the model as well and consider the accreting box model from Section 12.2.4 with delayed type Ia supernova enrichment. As we saw in that section, adding inflows of gas is a good way to solve the G dwarf problem and if inflows can also lead to better agreement on the \([\mathrm{O/Fe}]\) versus \([\mathrm{Fe/H}]\) distribution in the solar neighborhood, the inflow model would become even more compelling.
As in Section 12.2.4, we will assume that the gas mass is constant in time: \(M_g(t) = M_g(0)\), such that inflows balance the loss of gas due to outflows and the consumption of gas by star formation that is not returned through recycling. Because the gas mass is constant, the star formation rate is constant as well for a constant star formation efficiency time scale \(\tau_*\)
For a constant star formation rate, the DTD-averaged star formation rate follows from its definition in Equation \eqref{eq-chemev-DTDavgSFH}. For an exponentially declining DTD with decay timescale \(\tau_{\mathrm{Ia}}\) and minimum delay time \(\Delta t^{\mathrm{Ia}}_{\mathrm{min}}\)
The evolution equation of \(\dot{M}_{Z_\mathrm{Fe}^{\mathrm{Ia}}}\) is still given by Equation \eqref{eq-chemev-delayleaky-dmZfull-1} (we will immediately skip to the model that includes outflows and recycling) and plugging in the DTD-averaged star formation rate from Equation \eqref{eq-chemev-DTDMstar-accretingbox}, we get
or because \(\dot{M}_g = 0\)
We can solve this differential equation in the same way as we solved the spherical Jeans equation for arbitrary anisotropy profile \(\beta(r)\) in Chapter 6.4.2. We first compute the integrating factor
and then we can write Equation \eqref{eq-chemev-delayaccrete-dZfull-1} as
which we can directly integrate to give
where again \(\tilde{\tau}^{-1} = \tau_{\mathrm{Ia}}^{-1}-\tau_{\mathrm{dep}}^{-1}\) and \(t > \Delta t^{\mathrm{Ia}}_\mathrm{min}\).
To obtain the full evolution of iron and oxygen, we also need the equivalent expression for instantaneous recycling of type II supernovae enrichment. We could derive this starting from the basic evolution equation, but alternatively, we can simply take the limit of Equation \eqref{eq-chemev-ZFeIat-accretingbox} as \(\tau_{\mathrm{Ia}} \rightarrow 0\) and \(\Delta t^{\mathrm{Ia}}_\mathrm{min} \rightarrow 0\). In that limit, \(\tilde{\tau} \rightarrow \tau_{\mathrm{Ia}} \rightarrow 0\) as well and the expression therefore becomes simply, e.g., for the core-collapse contribution to the iron abundance
By working out the total mass as a function of time, it is straightforward to show that this expression is equivalent to that in Equation \eqref{eq-accretingbox-zM-outflow} (setting \(r=0\)): \(M_g\) is constant and \(M_* = M_g\,t/\tau_*\) from solving Equation \eqref{eq-chemev-dotmstar-delayedaccretingbox} with the boundary condition that \(M_* = 0\) at \(t=0\); we then have that \((1-M/M_g) = -M_*/M_g = -t/\tau_*\) and plugging this into Equation \eqref{eq-accretingbox-zM-outflow}, we recover Equation \eqref{eq-chemev-ZFeIIt-accretingbox}. Equation \eqref{eq-chemev-ZFeIIt-accretingbox} is also equivalent to Equation \eqref{eq-accretingbox-zM-outflow-recycling}, because in the presence of recycling, \((1-M/M_g) = -(1-r)\,t/\tau_*\).
To apply this model to the abundance distribution in the solar neighborhood, we need to set \(\tau_*\) and \(\eta\) (we’ll assume that \(r=0.4\), as derived from stellar evolution). From the observation that \(M/M_g = 10 \approx M_*/M_g = t/\tau_*\) after \(\approx 10\,\mathrm{Gyr}\) of evolution, we require that \(\tau_* \approx 1\,\mathrm{Gyr}\). This is a much smaller value than that required in the closed or leaky box models. But it is in better agreement with direct measurements of \(\tau_*\) in external galaxies, where both the gas mas \(M_g\) and the star formation rate \(\dot{M}_*\) can be directly determined observationally and ratio’d to obtain \(\tau_*\) (Leroy et al. 2008).
To set \(\eta\), we require that the oxygen abundance of the ISM is approximately solar, with \(Z_{\mathrm{O},\odot} = 0.0056\) from Equation \eqref{eq-chemev-solarO}. Using \(\tau_* = 1\,\mathrm{Gyr}\), we can then solve Equation \eqref{eq-chemev-ZFeIIt-accretingbox} applied to the oxygen abundance with \(p_\mathrm{O}^{\mathrm{II}} = 0.016\) for \(\eta\). The \(\eta\) parameter appears non-linearly, but the exponential term is essentially zero, so \(\eta \approx p_\mathrm{O}^{\mathrm{II}}/Z_{\mathrm{O},\odot}-1+r \approx 2.25\). Moderate outflows are therefore required to match the present day abundance of the ISM. Note that this derived value is higher the one we derived with similar reasoning in Section 12.2.4 for two reasons: (a) we now take into account recycling with \(r=0.4\); this increases the required \(\eta\) by \(r\); and (b) the ratio between the yield and the ISM abundance is higher when just applied to oxygen. We implement the accreting box model including outflows (so it is a leaky accreting box model) and compare it to the \([\mathrm{O/Fe}]\) and \([\mathrm{Fe/H}]\) data in the solar neighborhood:
[20]:
def feh_accretingbox(t,p_Fe_Ia=0.001,p_Fe_II=0.0011,
tau_Ia=0.4,tau_Ia_2=3.3,mindt=0.05,
tau_star=1.,solar_Fe=7.5,
eta=2.25,r=0.4):
tau_dep= tau_star/(1.+eta-r)
tau_tilde= 1./(1./tau_Ia-1./tau_dep)
dt= t-mindt
out= p_Fe_Ia/(1.+eta-r)\
*(1.-numpy.exp(-dt/tau_dep)
+tau_tilde/tau_dep*(numpy.exp(-dt/tau_Ia)
-numpy.exp(-dt/tau_dep)))
# Second exponential Ia component
if not tau_Ia_2 is None:
out/= 2.
tau_tilde= 1./(1./tau_Ia_2-1./tau_dep)
out+= p_Fe_Ia/2./(1.+eta-r)\
*(1.-numpy.exp(-dt/tau_dep)
+tau_tilde/tau_dep*(numpy.exp(-dt/tau_Ia_2)
-numpy.exp(-dt/tau_dep)))
out[t<mindt]= 0.
# Add CCSN contribution
out+= p_Fe_II/(1.+eta-r)*(1.-numpy.exp(-t/tau_dep))
# Convert to [Fe/H]
logZFe_solar= -2.90+solar_Fe-7.50
return numpy.log10(out)-logZFe_solar
def oh_accretingbox(t,p_O_II=0.016,
tau_star=1.,solar_O=8.69,
eta=2.25,r=0.4):
tau_dep= tau_star/(1.+eta-r)
# Just CCSN contribution
out= p_O_II/(1.+eta-r)*(1.-numpy.exp(-t/tau_dep))
# Convert to [O/H]
logZO_solar= -2.25+solar_O-8.69
return numpy.log10(out)-logZO_solar
ts= numpy.geomspace(1e-3,10.,201)
figsize(7,5)
plot_sn_ofe_feh(xrange=[-2.25,0.69],yrange=[-0.21,0.6])
plot(feh_closedbox(ts),oh_closedbox(ts)-feh_closedbox(ts),
lw=1.,label=r'$\mathrm{closed\ box}$')
plot(feh_accretingbox(ts),oh_accretingbox(ts)-feh_accretingbox(ts),
lw=2.,label=r'$\mathrm{accreting\ leaky\ box}$')
plot(feh_accretingbox(ts,p_Fe_Ia=0.0017,p_Fe_II=0.0012),
oh_accretingbox(ts)-feh_accretingbox(ts,p_Fe_Ia=0.0017,p_Fe_II=0.0012),
lw=2.,label=r'$\mathrm{accreting\ leaky\ box}; p_{\mathrm{Fe}}^{\mathrm{Ia}} = 0.0017\,, p_{\mathrm{Fe}}^{\mathrm{II}} = 0.0012$')
xlabel(r'$[\mathrm{Fe/H}]$')
ylabel(r'$[\mathrm{O/Fe}]$')
legend(frameon=False,fontsize=16);
We see that the model that we described above fits the higher \([\mathrm{O/Fe}]\) data well, quantitatively matching the behavior including the downturn at \([\mathrm{O/Fe}] \approx -0.7\). However, while the oxygen abundance of the ISM is approximately solar by design, the iron abundance of the ISM in the model at the present day is \(\approx -0.24\) and after 10 Gyr therefore \([\mathrm{O/Fe}] \approx 0.24\), while the data extend down to \([\mathrm{O/Fe}] \approx 0\). One reason for this may be the assumed yields, which are uncertain for iron, because (a) the type Ia yield depends on the normalization of the delay time distribution, which is uncertain at the \(\approx 50\%\) level (Maoz & Mannucci 2012) and the yield \(p_{\mathrm{Fe}}^{\mathrm{Ia}}\) could therefore by \(50\%\) higher leading to a higher present-day iron abundance; and (b) iron formed in core-collapse supernovae is produced during the explosive nucleosynthesis stage, whose outcome is more uncertain. The figure above therefore also shows a model with increased iron yields which get the present-day ISM abundance of iron up to \([\mathrm{Fe/H}] \approx -0.1\).
Another option is that gas is flowing in more slowly than in the specific accreting box model that has a constant gas mass that we considered here. If the gas mass decreases slowly over time (but not so slow as in the closed box model), then the present-day iron abundance of the ISM comes out higher, because the same yields from supernovae explosions get mixed into a smaller overall gas volume and thus increase the abundance more. In this case, the star-formation rate will be dropping over time rather than staying constant as in the model considered here. As discussed in Section 12.2.4, a slowly declining gas mass would also lead to better agreement with the observed metallicity distribution.
We therefore see that a galactic evolution model that includes realistic amounts of inflows and outflows of gas and the empirically determined delay time distribution of type Ia supernovae explosions, is able to approximately match both the metallicity distribution of stars (see Section 12.2.4) and the track of \([\mathrm{O/Fe}]\) versus \([\mathrm{Fe/H}]\) in the solar neighborhood. Uncertainties in the yields from stellar evolution (and the initial mass function) limit how well we can expect our models to agree with the data, certainly when considering more than these two relatively simple elements. And as is clear from the abundance distribution, the overall distribution has additional complexity in the form of two sequences of tracks in \([\mathrm{O/Fe}]\) versus \([\mathrm{Fe/H}]\) that are ubiquitous throughout the Milky Way (Hayden et al. 2015) and that cannot be captured by the single-zone models that we discussed in this chapter. But a simple star formation law, population-level stellar yields, inflows, outflows, and delays in type Ia supernova enrichment allows an overall understanding of galactic chemical evolution that is in good agreement with the observations.