5. Enrichment

VICE takes a general approach in modeling nucleosynthesis. All elements are treated equally; there are no special considerations for any element. In this documentation we derive the analytic form of the enrichment equation for an arbitrary element \(x\) with arbitrary nucleosynthetic yields for arbitrary evolutionary histories. This is an integro-differential equation of the element’s mass as a function of time, which VICE solves as an initial-value problem by imposing the boundary condition that its abundance at time zero is given by the primordial abundance from big bang nucleosynthesis. In this version of VICE, helium is the only element for which this value is nonzero.

5.1. The Enrichment Equation

The enrichment equation quantifies the rate of change of an element’s total mass present in the interstellar medium (ISM). At its core, it is a simple sum of source and sink terms.

\[\dot{M}_x = \dot{M}_x^\text{CC} + \dot{M}_x^\text{Ia} + \dot{M}_x^\text{AGB} - \frac{M_x}{M_g}\left[ \dot{M}_\star + \xi_\text{enh}\dot{M}_\text{out} \right] + \dot{M}_x^\text{r} + Z_{x,\text{in}}\dot{M}_\text{in}\]

where \(M_x\) is the mass of the element \(x\) in the interstellar medium, \(\dot{M}_x\) its time-derivative, and \(M_g\) the mass of the ISM gas. \(\dot{M}_x^\text{CC}\), \(\dot{M}_x^\text{Ia}\), and \(\dot{M}_x^\text{AGB}\) quantify the rate of production from core-collapse supernovae (CCSNe), type Ia supernovae (SNe Ia), and asymptotic giant branch (AGB) stars, respectively.

We detail each term individually here.

5.2. Core Collapse Supernovae

Core collapse supernovae (CCSNe) are the explosions of massive stars (\(\gtrsim 8 M_\odot\)) at the end of their post main sequence lifetimes. Due to the steep nature of the lifetime-stellar mass relationship, these stars have lifetimes that are extremely short compared to the relevant timescales of galactic chemical evolution (\(\sim\) few Myr compared to \(\sim\) few Gyr). To a good approximation, the lifetimes of these stars can be treated as instantaneous in zone models.

Note

Another motivation for this approximation is that the lifetimes are often significantly shorter than the typical mixing timescales in even modestly sized galaxies. The longest lifetimes of these stars is of order tens of megayears; in comparison, the mixing timescale in the solar annulus of the Milky Way is likely comparable to the dynamical timescale at this distance (\(\sim\) 250 Myr, a factor of ten larger). Zone models at their core already assume that these mixing timescales are negligibly short due to the assumption of instantaneous mixing; if CCSN timescales are even shorter, then they can certainly also be modeled as instantaneous.

VICE therefore approximates CCSNe as being simultaneous with the formation of their progenitor stars. This implies a linear relationship between the rate of production of some element \(x\) from CCSNe and the star formation rate:

\[\dot{M}_x^\text{CC} = \epsilon_x^\text{CC} y_x^\text{CC}(Z)\dot{M}_\star\]

where \(y_x^\text{CC}\) is the IMF-averaged fractional net yield of the element \(x\) from CCSNe at a metallicity \(Z\): the fraction of the entire stellar population’s initial mass that is processed into the element \(x\) and ejected to the interstellar medium minus the amount that the star was born with. \(\epsilon_x^\text{CC}\) is the entrainment fraction of the element \(x\) from CCSNe; this is the mass fraction of the net yield which is retained by the interstellar medium, the remainder of which is added directly to the outflow.

Note

VICE implements recycling of previously produced elements separate from nucleosynthesis, running from the standpoint of net rather than absolute yields.

In practice, \(y_x^\text{CC}\) is highly uncertain [1]. VICE therefore makes no assumptions about the user’s desired form of the yield; this parameter can be assigned either a number to represent a metallicity-independent yield, or a function of the metallicity by mass \(Z = M_x/M_g\). VICE includes features which will calculate the value of \(y_x^\text{CC}\) for a given element and metallicity based on the results of supernova nucleosynthesis studies upon request, but requires the user to specify an exact number or function.

Relevant Source Code:

  • vice/src/singlezone/ccsne.c

  • vice/core/dataframe/_yield_settings.pyx

  • vice/yields/ccsne/__init__.py

5.2.1. Extension to Multizone Models

Because VICE approximates core collapse supernovae as occuring instantaneously following the formation of their progenitor stars, this implies that CCSN progenitors also should not migrate between zones in a multizone model. Therefore, the formalism implemented for singlezone models is retained in multizone simulations.

Relevant Source Code:

  • vice/src/multizone/ccsne.c

5.3. Type Ia Supernovae

Type Ia supernovae are the thermonuclear detonations of white dwarf stars. Being the remnants of lower-mass stars, white dwarfs are born and explode on timescales longer than the mixing timescales of galaxies. Therefore, the intrinsic time delay is non-negigible.

This requires a model for the SN Ia delay-time distribution (DTD), defined as the rate of SN Ia explosions associated with a single stellar population. Given a DTD \(R_\text{Ia}\) and an age \(\tau\), the rate of production of some element \(x\) from a single stellar population is given by

\[\dot{M}_x^\text{Ia} = \epsilon_x^\text{Ia} y_x^\text{Ia}(Z)M_* \frac{ R_\text{Ia}(\tau) }{ \int_0^\infty R_\text{Ia}(t) dt }\]

Note

The integral of this equation from \(t = 0\) to \(\infty\) must equal the yield times the mass of the stellar population. This necessitates the normalization of the DTD.

where \(y_x^\text{Ia}\) is the IMF-averaged fractional net yield of the element \(x\) from SNe Ia at metallicity \(Z\): the fraction of the stellar population’s initial mass that is processed into the element \(x\) and ejected to the interstellar medium minus the amount that the star was born with. \(\epsilon_x^\text{Ia}\) is the entrainment fraction of the element \(x\) from SNe Ia; this is the mass fraction of the net yield which is retained by the interstellar medium, the remainder of which is added directly to the outflow.

Note

VICE implements recycling of previously produced elements separate from nucleosynthesis, running from the standpoint of net rather than absolute yields.

In practice, \(y_x^\text{Ia}\) is highly uncertain [2]. VICE therefore makes no assumptions about the user’s desired form of the yield; this parameter can be assigned either a number to represent a metallicity-independent yield or a function of metallicity by mass \(Z = M_x/M_g\). VICE includes features which will calculate the value of \(y_x^\text{Ia}\) for a given element and metallicity based on the results of supernova nucleosynthesis studies upon request, but requires the user to specify an exact number or function.

The rate of enrichment from all previous episodes of star formation can be derived by integrating this equation over all ages:

\[\dot{M}_x^\text{Ia} = \epsilon_x^\text{Ia} y_x^\text{Ia}(Z)\frac{ \int_0^t \dot{M}_*(t')R_\text{Ia}(t - t')dt' }{ \int_0^\infty R_\text{Ia}(t')dt' }\]

This can also be expressed as the star formation history up to a time \(t\) weighted by the SN Ia rate. VICE approximates this equation as:

\[\dot{M}_x^\text{Ia} \approx \frac{ \sum_i \epsilon_x^\text{Ia} y_x^\text{Ia}(Z_\text{ISM}(i\Delta t)) \dot{M}_*(i\Delta t) R_\text{Ia}(t - i\Delta t) \Delta t }{ \sum_i^{T_\text{Ia}} R_\text{Ia}(i\Delta t) \Delta t }\]

where the sum in the numerator is over all timesteps and in the denominator up to a time \(T_\text{Ia}\) denoting an adopted full length of the SN Ia duty cycle. The constant RIA_MAX_EVAL_TIME declares \(T_\text{Ia}\) = 15 Gyr in vice/src/sneia.h.

In implementation, VICE normalizes the DTD at the beginning of the simulation. For an age \(\tau = t - t'\):

\[R_\text{Ia}(\tau) \rightarrow \frac{ R_\text{Ia}(\tau) }{ \int_0^{T_\text{Ia}} R_\text{Ia}(\tau) d\tau } \approx \frac{ R_\text{Ia}(t - i\Delta t) }{ \sum_i^{T_\text{Ia}} R_\text{Ia}(i\Delta t)\Delta t } \implies R_\text{Ia}(t - t')\Delta t \rightarrow \frac{ R_\text{Ia}(t - i\Delta t)\Delta t }{ \sum_i^{T_\text{Ia}} R_\text{Ia}(i\Delta t)\Delta t }\]

Inserting the normalized rate into the equation for \(\dot{M}_x^\text{Ia}\):

\[\dot{M}_x^\text{Ia} \approx \sum_i \epsilon_x^\text{Ia} y_x^\text{Ia}(Z_\text{ISM}(i\Delta t)) \dot{M}_*(i\Delta t) R_\text{Ia}(t - i\Delta t)\]

VICE implements this normalization of \(R_\text{Ia}\) at the beginning of simulations due to the simplification of this expression introduced in doing so. This reduces the computational expense in calculating this quantity for each element at each timestep.

VICE includes two built-in DTDs, denoting by strings as plaw and exp. As their names suggest, they are a power-law and an exponential DTD:

  • “plaw”: \(R_\text{Ia} \sim t^{-1.1}\)

  • “exp”: \(R_\text{Ia} \sim e^{-t/\tau_\text{Ia}}\)

Users may also construct their own functional forms of \(R_\text{Ia}\), which must accept time in Gyr as the only parameter. These functions need not be normalized in any way; VICE normalizes the DTD automatically.

Relevant Source Code:

  • vice/src/sneia.h

  • vice/src/singlezone/sneia.c

  • vice/yields/sneia/__init__.py

5.3.1. Extension to Multizone Models

The migration of star particles into and out of zones can affect the SN Ia rate in a given zone. In a singlezone simulation it is exactly as expected for the star formation history, but in a multizone model, it is coupled to the star formation histories in other zones. Because VICE knows the zone each star particle occupies at all times in simulation, the rate of SN Ia production of some element \(x\) should be expressed not as an integral over the star formation history, but as a summation over the stellar populations present in the zone:

\[\dot{M}_x^\text{Ia} \approx \sum_i \epsilon_x^\text{Ia} y_x^\text{Ia}(Z_i) M_i R_\text{Ia}(\tau_i)\]

where \((Z_i)\), \(M_i\), and \(\tau_i\) are the metallicity, initial mass, and age, respectively, of the \(i\)’th star particle in a given zone at a given time.

It is important to note that with this implementation, VICE is not sampling the SN Ia delay time distribution and stochastically ejecting iron according thereto. Instead, star particles gradually eject iron to the interstellar medium with a time-dependence given by the delay-time distribution.

Relevant Source Code:

  • vice/src/multizone/sneia.c

5.4. Asymptotic Giant Branch Stars

Asymptotic giant branch (AGB) stars are evolved stars that have carbon-oxygen cores surrounded by helium and hydrogen shells. These stars undergo thermal pulsations due to explosive ignition of helium fusion in the shell, typically referred to as helium shell flashes. During these pulses, material from the core is often mixed into the outer layers via convection, a process known as dredge-up. This brings heavy nuclei produced in the deeper regions of the star to the envelope, which is then ejected to the interstellar medium (ISM). This is one of the primary sites of s-process nucleosynthesis in the universe.

It may be tempting to model AGB star enrichment as a delay-time distribution (DTD) similar to that adopted for SNe Ia. However, this approach would implicitly adopt the assumption that every element is enriched via AGB stars with the same DTD, or that for a given element, the effective DTD is independent of metallicity. These may be fine assumptions, but it is not adopted in VICE due to the desire for as few assumptions as possible.

Instead, AGB star enrichment in VICE is implemented using the mass-lifetime relationship for stars and the main sequence mass fraction (MSMF). However, the form of the MSMF required here differs in detail from the true MSMF. Being evolved stars, the MSMF does not consider AGB stars. It is thus not the MSMF and the main sequence lifetimes of stars that are of interest, but the mass fraction of both main sequence and evolved stars and the total lifetime of stars. The form of \(h(t)\) necessary for modeling AGB star enrichment then changes to:

\[h(t) \rightarrow \frac{ \int_l^{m_\text{postMS}(t)} m \frac{dN}{dm} dm }{ \int_l^u m \frac{dN}{dm} dm }\]

The numerator is evaluated from \(l\) to the mass of stars ending their post main sequence lifetime \(m_\text{postMS}\) rather than the main sequence turnoff mass \(m_\text{to}\). As detailed here for a stellar population of age \(\tau\):

\[m_\text{postMS} = \left(\frac{\tau}{(1 + p_\text{MS})\tau_\odot} \right)^{-1/\alpha}\]

where \(\alpha\) is the power-law index on the mass-lifetime relationship, \(\tau_\odot\) is the main sequence lifetime of the sun, and \(p_\text{MS}\) is the ratio of a star’s post main sequence lifetime to its main sequence lifetime.

From a single stellar population, the rate of ejection of an element \(x\) from AGB stars to the ISM is given by:

\[\dot{M}_x^\text{AGB} = -\epsilon_x^\text{AGB}y_x^\text{AGB}(m_\text{postMS}, Z)M_\star\dot{h}\]

where \(\dot{h}\) is evaluated at the lookback time to the stellar population’s formation [3], \(M_\star\) is the initial mass of the stellar population, and \(y_x^\text{AGB}\) is the fractional net yield of \(x\) from an AGB star of initial mass \(m_\text{postMS}\) and metallicity \(Z\): the fraction of a single star’s initial mass that is processed into element \(x\) and ejected to the interstellar medium minus the amount that the star was born with. \(\epsilon_x^\text{AGB}\) is the entrainment fraction of the element \(x\) from AGB stars; this is the mass fraction of the net yield which is retained by the interstellar medium, the remainder of which is added directly to the outflow.

Note

VICE implements recycling of previously produced elements separate from nucleosynthetic yields, running from the standpoint of net rather than absolute yields.

For continuous star formation, the enrichment rate can be expressed as this quantity integrated over the star formation history:

\[\dot{M}_x^\text{AGB} = -\int_0^t \epsilon_x^\text{AGB} y_x^\text{AGB}(m_\text{postMS}(t - t'), Z_\text{ISM}(t')) \dot{M}_\star(t') \dot{h}(t - t') dt\]

This expression is approximated numerically as:

\[\dot{M}_x^\text{AGB} \approx \sum_i \epsilon_x^\text{AGB} y_x^\text{AGB}(m_\text{postMS}(t - i\Delta t), Z_\text{ISM}(i\Delta t)) \dot{M}_\star(i\Delta t) \left[h((i + 1)\Delta t) - h(i\Delta t)\right]\]

where the summation is taken over all previous timesteps. The need to differentiate \(h\) with time is eliminated in the numerical approximation by allowing each stellar population to be weighted by \(\Delta h\) between the current timestep and the next, made possible by the quantization of timesteps.

In practice, \(y_x^\text{AGB}\) is highly uncertain [4]. VICE therefore makes no assumptions about the user’s desired form of the yield; this parameter can be assigned either a built-in table published in an AGB star nucleosynthesis study or a function of stellar mass and metallicity constructed by the user.

Relevant source code:

  • vice/src/singlezone/agb.c

  • vice/core/dataframe/_agb_yield_settings.pyx

  • vice/yields/agb/__init__.py

5.4.1. Extension to Multizone Models

The migration of star particles into and out of zones can affect the AGB star enrichment rate in a given zone. In a singlezone simulation it is exactly as expected for the star formation history, but in a multizone model, it is coupled to the star formation histories in other zones. Because VICE knows the zone each star particle occupies at all times in simulation, the rate of AGB star enrichment rate of some element \(x\) should not be expressed as an integral over the star formation history, but as a summation over the stellar populations present in the zone:

\[\dot{M}_x^\text{AGB} \approx \sum_i \epsilon_x^\text{AGB} y_x^\text{AGB}(m_\text{postMS}(\tau_i), Z_i) M_i [h(\tau_i) - h(\tau_i + \Delta t)]\]

where \(Z_i\), \(M_i\), and \(\tau_i\) are the metallicity, initial mass, and age, respectively, of the \(i\)’th star particle in a given zone at a given time.

Relevant Source Code:

  • vice/src/multizone/agb.c

5.5. Subsequent Terms

The remaining terms in the enrichment equation make simple statements about remaining source and sink terms.

VICE retains the assumption that stars are born at the same metallicity as the ISM from which they form. This motivates the sink term

\[-\left(\frac{M_x}{M_g}\right)\dot{M}_\star\]

where the mass of the element \(x\) is depleted at the metallicity of the ISM \(Z_x = M_x/M_g\) in proportion with the star formation rate \(\dot{M}_\star\).

Many galactic chemical evolution models to date have assumed that outflows from galaxies occur at the same metallicity of the ISM. This would suggest that \(\dot{M}_x^\text{out} \approx (M_x/M_g)\dot{M}_\text{out}\). However, recent work in the astronomical literature from both simulations (e.g. Christensen et al. (2018) [5]) and observations (e.g. Chisholm, Trimonti & Leitherer (2018) [6]) suggest that this may not be the case. Therefore, VICE allows outflows to occur at some multiplicative factor \(\xi_\text{enh}\) above or below the ISM metallicity, which may vary with time. This motivates the sink term

\[-\left(\frac{M_x}{M_g}\right)\xi_\text{enh}\dot{M}_\text{out}\]

Because VICE works with net rather than absolute yields, simulations must quantify the rate at which stars return mass to the ISM at their birth metallicity. This is mathematically similar to the rate of total gas recycling, but weighted by the metallicities of the stars recycling. Since stars are assumed to form at the metallicity of the ISM,

\[\dot{M}_x^\text{r} = \int_0^t \dot{M}_\star(t') Z_{x,\text{ISM}}(t') \dot{r}(t - t') dt\]

where \(r(\tau)\) is the cumulative return fraction from a single stellar population of age \(\tau\). This is approximated numerically as

\[\dot{M}_x^\text{r} \approx \sum_i \dot{M}_\star(i\Delta t) Z_{x,\text{ISM}}(i\Delta t) \left[r((i + 1)\Delta t) - r(i\Delta t)\right]\]

where the summation is taken over all previous timesteps. The need to differentiate \(r\) with time is eliminated in the numerical approximation by allowing each stellar population to be weighted by \(\Delta r\) between the current timestep and the next, made possible by the quantization of timesteps. In the event that the user has specified instantaneous recycling:

\[\dot{M}_x^\text{r} = r_\text{inst}\dot{M}_\star Z_{x,\text{ISM}}\]

At any given timestep, there is gas infall onto the simulated galaxy of a given metallicity \(Z\). In most cases this term is negligibly small, but in some interesting cases it may not be (e.g. a major merger event). This necessitates the final term \(Z_{x,\text{in}}\dot{M}_\text{in}\).

Relevant Source Code:

  • vice/src/singlezone/recycling.c

  • vice/src/singlezone/element.c

  • vice/src/singlezone/ism.c

5.5.1. Extension to Multizone Models

The only subsequent term of the enrichment equation modified in multizone simulations is that quantifying the rate of recycling of an element \(x\). The migration of star particles into and out of zones can affect the recycling rate in a given zone. In a singlezone simulation it is exactly as expected for the star formation history, but in a multizone model, it is coupled to the star formation histories in other zones. Because VICE knows the zone each star particle occupies at all times in simulation, the rate of recycling of some element \(x\) should be expressed not as an integral over the star formation history, but as a summation over the stellar populations in the zone:

\[\dot{M}_x^\text{r} \approx \sum_i M_i Z_{x,i} [r(\tau_i + \Delta t) - r(\tau_i)]\]

where \(M_i\), \(Z_i\), and \(\tau_i\) are the mass, metallicity, and age, respectively, of the \(i\)’th star particle in a given zone at a given time.

Relevant Source Code:

  • vice/src/multizone/recycling.c

  • vice/src/multizone/element.c

  • vice/src/multizone/ism.c

5.6. Sanity Checks

At all timesteps VICE forces the mass of every element to be non-negative. If the mass is found to be below zero at any given time, it is assumed to not be present in the interstellar medium and is assigned a mass of exactly zero. Absent this, the mass of each element reported by VICE is merely the numerically estimated solution to the enrichment equation.

Relevant source code:

  • vice/src/singlezone/element.c