On the deuterium abundance and the importance of stellar mass loss in the interstellar and intergalactic medium

TL;DR

This study uses cosmological simulations to analyze deuterium abundance, revealing its close relation to primordial levels, the impact of stellar mass loss, and providing predictions for future observations across different galactic environments.

Contribution

It offers a detailed quantification of deuterium evolution in galaxies, emphasizing the role of stellar mass loss and comparing simulation results with high-redshift observations.

Findings

Deuterium abundance remains near primordial at low metallicity.

Deuterium fraction decreases slightly with redshift due to stellar mass loss.

Simulations match observed deuterium levels in the Milky Way's interstellar medium.

Abstract

We quantify the gas-phase abundance of deuterium and fractional contribution of stellar mass loss to the gas in cosmological zoom-in simulations from the Feedback In Realistic Environments project. At low metallicity, our simulations confirm that the deuterium abundance is very close to the primordial value. The chemical evolution of the deuterium abundance that we derive here agrees quantitatively with analytical chemical evolution models. We furthermore find that the relation between the deuterium and oxygen abundance exhibits very little scatter. We compare our simulations to existing high-redshift observations in order to determine a primordial deuterium fraction of 2.549 +/- 0.033 x 10^-5 and stress that future observations at higher metallicity can also be used to constrain this value. At fixed metallicity, the deuterium fraction decreases slightly with decreasing redshift, due to the increased importance of mass loss from intermediate-mass stars. We find that the evolution of the average deuterium fraction in a galaxy correlates with its star formation history. Our simulations are consistent with observations of the Milky Way's interstellar medium: the deuterium fraction at the solar circle is 85-92 per cent of the primordial deuterium fraction. We use our simulations to make predictions for future observations. In particular, the deuterium abundance is lower at smaller galactocentric radii and in higher mass galaxies, showing that stellar mass loss is more important for fuelling star formation in these regimes (and can even dominate). Gas accreting onto galaxies has a deuterium fraction above that of the galaxies' interstellar medium, but below the primordial fraction, because it is a mix of gas accreting from the intergalactic medium and gas previously ejected or stripped from galaxies.

Equations

$$f_j=\frac{\frac{m_j}{{\rho }_j}W(r_j,h_{\mathrm{sml}})}{{\Sigma }_i\frac{m_i}{{\rho }_i}W(r_i,h_{\mathrm{sml}})},$$

$$\begin{array}{rl}\frac{(D/H)}{{(D/H)}_{\mathrm{prim}}}& =\frac{\mathrm{D}}{{\mathrm{D}}_{\mathrm{prim}}}\frac{{\mathrm{H}}_{\mathrm{prim}}}{\mathrm{H}}=\frac{m_{\mathrm{initial}}}{m_{\mathrm{initial}}+m_{\mathrm{recycled}}}\frac{X_{\mathrm{prim}}}{X_{\mathrm{gas}}}\\ & =\frac{m_{\mathrm{initial}}}{m_{\mathrm{gas}}}\frac{X_{\mathrm{prim}}}{X_{\mathrm{gas}}},\end{array}$$

$$\begin{array}{rl}{f}_{\mathrm{recycled}}& =\frac{m_{\mathrm{recycled}}}{m_{\mathrm{gas}}}=\frac{m_{\mathrm{gas}}-m_{\mathrm{initial}}}{m_{\mathrm{gas}}}\\ & =1-\frac{\mathrm{D}}{{\mathrm{D}}_{\mathrm{prim}}}=1-\frac{(D/H)}{{(D/H)}_{\mathrm{prim}}}\frac{X_{\mathrm{gas}}}{X_{\mathrm{prim}}},\end{array}$$

$$\frac{\mathrm{D}}{{\mathrm{D}}_{\mathrm{prim}}}=\frac{1}{1+r{Z}_{\mathrm{O}}/m_{\mathrm{O}}},$$

$$\frac{\mathrm{D}}{{\mathrm{D}}_{\mathrm{prim}}}=e^{\frac{-r{Z}_{\mathrm{O}}}{m_{\mathrm{O}}}}$$

Full text

††pagerange: On the deuterium abundance and the importance of stellar mass loss in the interstellar and intergalactic medium–References††pubyear: 2018

1 Introduction

Deuterium is one of the few stable isotopes produced in astrophysically interesting amounts during Big Bang nucleosynthesis, together with helium and lithium (see Steigman 2007 for a review). Helium and lithium can be produced after this initial phase, in stars and via collisions of cosmic ray nuclei, potentially increasing their gas-phase abundances. However, the gas-phase deuterium abundance can only decrease. All primordial deuterium is burned during the collapse of a protostar and deuterium synthesized in stellar interiors is immediately destroyed, because deuterium fuses at relatively low temperatures, $T\approx 10^{6}\leavevmode\nobreak\ K$, easily reached in the interiors of stars and even brown dwarfs (Epstein et al., 1976; Stahler, 1988; Spiegel et al., 2011). Therefore, mass lost from stars (also referred to as ‘recycled’ gas) is deuterium-free, i.e. $\mathrm{(D/H)}_{\mathrm{recycled}}=0$.

The primordial deuterium fraction, $\mathrm{(D/H)}_{\mathrm{prim}}$, is sensitive to cosmological parameters and, in particular, to the baryon–photon ratio and thus to the baryonic density of the Universe. Measurements of the cosmic microwave background (CMB) radiation have pinned down the ratio of the mean density of baryons to the critical density and the Hubble parameter (Planck Collaboration XIII, 2016). The most recent theoretical models of Big Bang nucleosynthesis have incorporated these and derived, for example, $\mathrm{(D/H)}_{\mathrm{prim}}=(2.45\pm 0.05)\times 10^{-5}$ (Coc et al., 2015) and $\mathrm{(D/H)}_{\mathrm{prim}}=(2.58\pm 0.13)\times 10^{-5}$ (Cyburt et al., 2016), where the quoted errors are $1\sigma$.

An accurate determination of the primordial deuterium fraction, in conjunction with Big Bang nucleosynthesis reaction rates, gives an independent constraint on the cosmic baryon density. If there is disagreement with the value derived from CMB measurements, this could point to a deviation in the expansion rate of the early universe and to non-standard models of big bang nucleosynthesis. Low-metallicity gas likely has a deuterium fraction close to the primordial value, because it has not been substantially enriched by stellar mass loss. Absorption lines in spectra of background quasars have been used to determine the primordial deuterium fraction observationally, finding e.g. $\mathrm{(D/H)}_{\mathrm{prim}}=(2.547\pm 0.033)\times 10^{-5}$ (Cooke et al., 2016). Modern estimates are thus consistent with each other and there is currently no conflict with the standard model of cosmology (e.g. Steigman, 2007; Cooke et al., 2014).

Intermediate-mass and massive stars return material to the interstellar medium (ISM) via stellar winds before and during the asymptotic giant branch (AGB) phase and via supernova explosions, respectively. One well-known effect of this recycling process of baryons that become part of a star and are later returned into space is the release of metals into the ISM and the intergalactic medium (IGM). However, it is also important for the destruction of light elements, such as deuterium. If there is no fresh infall of gas onto galaxies and the ISM of these objects is replenished by stellar mass loss, both the metallicity of the gas and young stars increases and the deuterium fraction in the ISM decreases. The ratio of the deuterium fraction in the ISM or IGM and the primordial deuterium fraction, $\mathrm{(D/H)}/\mathrm{(D/H)}_{\mathrm{prim}}$, is therefore a measure of the fraction of the gas that has not been processed in stars. The inverse of this, i.e. $\mathrm{(D/H)}_{\mathrm{prim}}/\mathrm{(D/H)}$, is known as the astration factor. Measurements of the evolution of the deuterium abundance have been used to constrain galactic chemical evolution models (e.g. Audouze & Tinsley, 1974; Vangioni-Flam & Audouze, 1988; Vangioni-Flam et al., 1994; Scully et al., 1997; Olive et al., 2012). These models predict astration factors higher than observed in the local ISM when they only take into account cosmological inflow (e.g. Fields, 1996; Romano et al., 2006; Lagarde et al., 2012). Models that additionally allow for galactic outflows predict lower astration factors (e.g. Dvorkin et al., 2016; Weinberg, 2017). In this way, measurements of the deuterium fraction (and thus the astration factor) can shed light on the balance between primordial inflow, metal-enriched outflow, and recycling through stellar mass loss, which are all related to the star formation and accretion history of a galaxy (e.g. Cassé et al., 1998; Prantzos & Ishimaru, 2001; Romano et al., 2006; Dvorkin et al., 2016; Weinberg, 2017).

The fuelling of star formation by stellar mass loss is likely more important in high-mass galaxies and high-density environments. Massive early-type galaxies and satellite galaxies have specific star formation rates (SFRs) far below those of central late-type galaxies. It is not known which process(es) quench(es) galaxies, but galactic outflows are at least partially responsible for quenching massive galaxies and preventing subsequent gas accretion (e.g. Faucher-Giguère et al., 2011; van de Voort et al., 2011). However, a substantial fraction of local early-type galaxies still have a detectable molecular or atomic gas reservoir (e.g. Young et al., 2011; Serra et al., 2012). Some of these exhibit gas kinematics indicating a predominantly external gas supply, such as through minor mergers, whereas others (especially those located in a cluster environment) are consistent with their ISM being fed through stellar mass loss (Davis et al., 2011). Furthermore, some massive galaxies in the centres of clusters are forming stars at a substantial rate ($1-100$ M*☉* yr*-1*) and contain a considerable amount of dust (O’Dea et al., 2008; Donahue et al., 2011). Dust is produced by stars and destroyed by sputtering in hot gas. Therefore, the gas supply is unlikely to have cooled out of the hot halo gas. This also indicates that stellar mass loss may be an important contributor to the fuel for the observed star formation (Voit & Donahue, 2011).

The amount of mass supplied to the ISM through stellar mass loss could also be sufficient to fuel most of the star formation in present-day star-forming galaxies, including the Milky Way (Leitner & Kravtsov, 2011). However, galactic outflows were not included, but are likely required to produce correct stellar masses and metallicities. Hydrodynamical simulations that included feedback from stars and/or black holes have found that stellar mass loss becomes more important for fuelling star formation towards lower redshift, although, in general, it does not become the dominant fuel source for star formation (Oppenheimer & Davé, 2008; Segers et al., 2016).

The predicted deuterium fraction and the importance of stellar mass loss are the focus of this paper. We present results from a suite of high-resolution, cosmological ‘zoom-in’ simulations from the ‘Feedback In Realistic Environments’ (FIRE) project,111http://fire.northwestern.edu/ which spans a large range in halo and galaxy mass. The FIRE simulation suite has been shown to successfully reproduce a variety of observations, which is linked to the strong stellar feedback implemented. These galactic winds efficiently redistribute gas from galaxies out to large galactocentric distances (see Muratov et al. 2015, 2017). For the purposes of this paper, we highlight the fact that the simulations match the derived stellar–to–halo mass relationship (Hopkins et al., 2014; Feldmann et al., 2016), the galaxy mass-metallicity relation and gas-phase metallicity gradients at $z=0-3$ (Ma et al., 2016, 2017), and the dense neutral hydrogen, H i, content of galaxy haloes (Faucher-Giguère et al., 2015, 2016; Hafen et al., 2017).

This is the first time cosmological, hydrodynamical simulations are used to study the deuterium abundance in the ISM and IGM. Our simulations self-consistently follow the time-dependent assembly of dark matter haloes, the accretion of gas onto galaxies, the formation of stars, the return of mass in stellar winds, and the generation of large-scale galactic outflows, whereas chemical evolution models - which are often used for similar studies - are based on analytic prescriptions. Specifically, we do not assume instantaneous recycling of stellar mass loss nor instantaneous mixing of metals nor a specific parametrized gas accretion or star formation history. However, our simulations do not consider mixing of the gas between resolution elements.

In Section 2 we describe the suite of simulations used, as well as the way we compute the deuterium abundance and the fractional contribution of stellar mass loss to the gas, i.e. the ‘recycled gas fraction’ or $f_{\mathrm{recycled}}$ (Section 2.1). The deuterium retention fraction and recycled gas fraction are related via $\mathrm{(D/D_{prim})}=1-f_{\mathrm{recycled}}$. In Section 3 we present our results, including comparisons to existing observations. Section 3.1 describes the evolution of the deuterium fraction (and hence of the recycled gas fraction), while Section 3.2 focuses on high redshift and Section 3.3 on low redshift. We discuss our results and conclude in Section 4.

2 Method

The simulations used are part of the FIRE-1 sample. These were run with gizmo222http://www.tapir.caltech.edu/$\scriptstyle\sim$phopkins/Site/GIZMO.html (Hopkins, 2015) in ‘P-SPH’ mode, which adopts the Lagrangian ‘pressure-energy’ formulation of the smoothed particle hydrodynamics (SPH) equations (Hopkins, 2013). The gravity solver is a heavily modified version of gadget-2 (Springel, 2005), with adaptive gravitational softening following Price & Monaghan (2007). Our implementation of P-SPH also includes substantial improvements in the artificial viscosity, entropy diffusion, adaptive timestepping, smoothing kernel, and gravitational softening algorithm.

The FIRE project consists of a suite of cosmological ‘zoom-in’ simulations of galaxies with a wide range of masses, simulated to $z=0$ (Hopkins et al. 2014; Chan et al. 2015; Ma et al. 2016; Hafen et al. 2017; Feldmann et al. in preparation), to $z=1.7$ (Feldmann et al., 2016), and to $z=2$ (Faucher-Giguère et al., 2015). The simulation sample used is identical to the one used in van de Voort et al. (2016) and the simulation details are fully described in Hopkins et al. (2014) and references therein. The three Milky Way-mass galaxies that are the focus of Figure 3 and 5 are simulations ‘m12i’, ‘m12v’, and ‘m11.9a’ (from highest to lowest stellar mass) from Hopkins et al. (2014) and Hafen et al. (2017). A $\Lambda$CDM cosmology is assumed with parameters consistent with the 9-yr Wilkinson Microwave Anisotropy Probe (WMAP) results (Hinshaw et al., 2013). The initial particle masses for baryons (dark matter) vary from $2.6\times 10^{2}-4.5\times 10^{5}$ M*☉* ($1.3\times 10^{3}-2.3\times 10^{6}$ M*☉) for the 16 simulations that were run to $z=0$ (see also van de Voort et al. 2016 for further details). The 23 simulations that were run to $z\approx 2$ are described in Faucher-Giguère et al. (2015) and Feldmann et al. (2016) and their initial baryonic (dark matter) masses are $(3.3-5.9)\times 10^{4}$ M☉* ($(1.7-2.9)\times 10^{5}$ M*☉*).

Star formation is restricted to molecular, self-gravitating gas above a hydrogen number density of $n_{\mathrm{H}}\approx 5-50$ cm*-3*, where the molecular fraction is calculated following Krumholz & Gnedin (2011) and the self-gravitating criterion following Hopkins et al. (2013). The majority of stars form at gas densities significantly higher than this imposed threshold. Stars are formed from gas satisfying these criteria at the rate $\dot{\rho}_{\star}=\rho_{\mathrm{molecular}}/t_{\mathrm{ff}}$, where $t_{\mathrm{ff}}$ is the free-fall time. When selected to undergo star formation, the entire gas particle is converted into a star particle.

We obtain stellar evolution results from STARBURST99 (Leitherer et al., 1999) and assume an initial stellar mass function (IMF) from Kroupa (2002). Radiative cooling and heating are computed in the presence of the CMB radiation and the ultraviolet (UV)/X-ray background from Faucher-Giguère et al. (2009). Self-shielding is accounted for with a local Sobolev/Jeans length approximation. We impose a temperature floor of 10 K or the CMB temperature.

The primordial abundances are $X_{\mathrm{prim}}=0.76$ and $Y_{\mathrm{prim}}=0.24$, where $X_{\mathrm{prim}}$ and $Y_{\mathrm{prim}}$ are the mass fractions of hydrogen and helium, respectively. The simulations include a metallicity floor at metal mass fraction $Z_{\mathrm{prim}}\approx 10^{-4}$ Z*☉* or $Z_{\mathrm{prim}}\approx 10^{-3}$ Z*☉*, because yields are very uncertain at lower metallicities and we do not resolve the formation of individual first-generation stars. The abundances of 11 elements (H, He, C, N, O, Ne, Mg, Si, S, Ca and Fe) produced by massive and intermediate-mass stars are computed following Iwamoto et al. (1999), Woosley & Weaver (1995), and Izzard et al. (2004). The amount of mass and metals ejected in a computational time-step depends on the age of the star particle and our simulations therefore self-consistently follow time-dependent chemical enrichment. Mass ejected through supernovae and stellar winds are modelled by transferring a fraction of the mass of a star particle to its neighbouring gas particles, $j$, within its SPH smoothing kernel as follows:

[TABLE]

where $h_{\mathrm{sml}}$ is the smoothing length of the star particle (determined in the same manner as for gas particles), $r_{i}$ is the distance from the star particle to neighbour $i$, $W$ is the quintic SPH kernel, and the summation is over all SPH neighbours of the star particle, 62 on average. There is no sub-resolution metal diffusion in these simulations.

The FIRE simulations include an explicit implementation of stellar feedback by supernovae, radiation pressure, stellar winds, and photo-ionization and photo-electric heating (see Hopkins et al. 2014 and references therein for details). Feedback from active galactic nuclei (AGN) is not included. For star-forming galaxies, which constitute the majority of our simulated galaxies, AGN are thought to be unimportant. However, AGN-driven outflows are potentially important for the high-mass end of our simulated mass range.

We measure a galaxy’s stellar mass, $M_{\mathrm{star}}$, within 20 proper kpc of its centre. The deuterium fraction of a galaxy’s ISM is measured within 20 proper kpc of its centre for gas with a temperature below $10^{4}$ K, which selects the warm ionized and cold neutral gas in the ISM. These choices have a mild effect on the normalization of some of our results, but not on the trends or on our conclusions.

2.1 Deuterium fraction in hydrodynamical simulations

Determining $\mathrm{(D/H)}/\mathrm{(D/H)}_{\mathrm{prim}}$ in our simulations is straightforward. The mass of a gas particle can only increase during the simulation by receiving mass lost from nearby stars (no particle splitting is implemented). Therefore, any mass above the initial particle mass, $m_{\mathrm{initial}}$, is deuterium-free. This is mixed with the initial particle mass, which has the primordial fraction of deuterium. Therefore, for each gas particle, we calculate

[TABLE]

where $m_{\mathrm{recycled}}$ is the amount of mass received from evolving stars, i.e. the amount of gas that has been ‘recycled’, $X_{\mathrm{gas}}$ is the mass fraction of hydrogen, and $m_{\mathrm{gas}}$ is the mass of the particle at the redshift of interest. We refer to this quantity as the deuterium retention fraction, because it is the fraction of deuterium, produced during Big Bang nucleosynthesis, that is not destroyed. The inverse of Equation 2 is the astration factor. The value of $\mathrm{(D/H)}_{\mathrm{prim}}$ is well-constrained, both directly from absorption-line observations of low-metallicity gas and indirectly from CMB measurements coupled with Big Bang nucleosynthesis reaction rates (e.g. Cooke et al., 2016; Coc et al., 2015; Cyburt et al., 2016). Another way to constrain the primordial value is by comparing observations to cosmological simulations, as done in Section 3.2.

The deuterium retention fraction in Equation 2 is directly related to the fractional contribution of stellar mass loss to the gas, i.e. the recycled gas fraction,

[TABLE]

which is used to study the importance of stellar mass loss in fuelling the ISM and star formation. Besides destroying all deuterium, a fraction of the hydrogen is fused into helium and metals before the gas is recycled into the ISM, $X_{\mathrm{gas}}=1-Y_{\mathrm{gas}}-Z_{\mathrm{gas}}$. This can be approximated well by $X_{\mathrm{gas}}\approx X_{\mathrm{prim}}-3Z_{\mathrm{gas}}$. The factor $X_{\mathrm{gas}}/X_{\mathrm{prim}}$ is close to unity for subsolar metallicities, but becomes more important at supersolar metallicities. Even though the differences do not change our conclusions, we will show and discuss both $\mathrm{(D/H)}/\mathrm{(D/H)}_{\mathrm{prim}}$ and $\mathrm{D/D_{prim}}$ or $f_{\mathrm{recycled}}$ when relevant.

3 Results

Observations of the deuterium fraction exist at both high and low redshift. We will first discuss the evolution of (D/H) and then discuss predictions and observational comparisons at $z=3$ and $z=0$ separately. Throughout the paper, we use oxygen abundance ratios of gas as compared to those of the Sun, i.e. $\mathrm{[O/H]}=\mathrm{log_{10}}(n_{\mathrm{O}}/n_{\mathrm{H}})-\mathrm{log_{10}}(n_{\mathrm{O}}/n_{\mathrm{H}})_{\astrosun}$, where $n_{\mathrm{O}}$ is the oxygen number density, $n_{\mathrm{H}}$ the hydrogen number density, and $\mathrm{log_{10}}(n_{\mathrm{O}}/n_{\mathrm{H}})_{\astrosun}=-3.31$ is the solar oxygen abundance taken from Asplund et al. (2009).

3.1 Evolution of the deuterium fraction

While the total deuterium content of the Universe decreases with time, its total metallicity increases, leading to an inverse correlation between the deuterium and oxygen abundance (e.g. Steigman, 2003; Romano et al., 2006; Dvorkin et al., 2016). Figure 1 shows the median deuterium retention fraction (black curves) for all gas particles in our simulations as a function of oxygen metallicity at $z=0$ (top panels) and $z=3$ (bottom panels). The left panels only take into account deuterium and oxygen, whereas the right panels fold in the hydrogen abundance for both axes. $Z_{\mathrm{O}}$ ($Z_{\mathrm{O,\astrosun}}$) is the oxygen mass fraction of the gas (for solar abundances). The grey, shaded regions show the $1\sigma$ (dark) and $3\sigma$ (light) scatter around the median. Some of our simulations implemented a relatively high metallicity floor of $\mathrm{[O/H]}_{\mathrm{initial}}=-2.8$. Here, in order to not be affected by the imposed metallicity floor, we excluded gas with a metallicity within a factor of 2 from its initial oxygen abundance, but this choice does not affect our conclusions. At solar oxygen metallicity, about 90 per cent of the primordial deuterium is not destroyed at $z=3$ and 88 per cent at $z=0$. $\mathrm{(D/H)}/\mathrm{(D/H)}_{\mathrm{prim}}=0.91$ ($z=0$) and $0.93$ ($z=3$) at $\mathrm{[O/H]}=0$, slightly higher because of the small decrease of the hydrogen fraction.

The $1\sigma$ scatter in this relation is very small, which shows that the destruction of deuterium and the enrichment with oxygen are tightly correlated. However, the scatter increases at $z=0$ at the highest metallicities ($\mathrm{[O/H]}\gtrsim 0.5$). Additionally, we find large non-Gaussian tails at all metallicities, which means that even at low metallicity, a small fraction of gas particles have substantially reduced deuterium abundances. Our calculations likely underestimate the mixing of gas, because elements in our simulation are stuck to gas particles and do not diffuse to neighbouring gas particles. Adding turbulent diffusion to our simulations would only decrease the scatter in the correlation between deuterium and oxygen, because it smoothes out variations, and would thus strengthen our conclusions. The dependence of (D/H) on [O/H] is very steep at high metallicity, because [O/H] is a logarithmic quantity.

A small fraction of the gas (0.5 per cent) reaches extremely high metallicities ($\mathrm{[O/H]}>0.5$), which have not been observed. It is possible that such rare systems exist, outside the Milky Way, but are beyond current observational capabilities. Note, however, that the (average) metallicity in sightlines through our simulation is always $\mathrm{[O/H]}\leq 0.5$ (see Figure 4). Another possibility is that there is not enough mixing in our simulations, since the metals are stuck to particles and cannot diffuse, resulting in small metal-rich pockets. Additionally, the yields are very uncertain at such high metallicities. The real uncertainty is therefore larger than the scatter in this regime. Note that although $\mathrm{D/D_{prim}}\leq 1$, $\mathrm{(D/H)}/\mathrm{(D/H)}_{\mathrm{prim}}$ can be larger than unity in rare cases at high metallicity, because $\mathrm{H/H_{prim}}$ can become very small due to hydrogen fusion.

For comparison, the red, dashed curves (identical in top and bottom left panels) show the relation between the oxygen and deuterium abundances obtained from a one-zone chemical evolution model (Weinberg, 2017). This model assumes that chemical equilibrium is reached in the ISM due to the balance between gas inflow and outflow, enrichment though stellar mass loss, and gas consumption due to star formation. The only parameters in the relation are the recycling fraction, $r$, i.e. the fraction of mass returned to the ISM by a simple stellar population, and the oxygen yield, $m_{\mathrm{O}}$, i.e. the mass fraction of a simple stellar population released into the ISM in oxygen,

[TABLE]

where $Z_{\mathrm{O}}$ is the oxygen mass fraction of the gas. We also compare our findings to analytic results derived from a closed box model, i.e. no gas inflow or outflow (Tinsley, 1980). The resulting relation between the deuterium retention fraction and oxygen abundance is333Note that in the literature, this equation is usually expressed with a yield defined with respect to the final stellar remnant mass (after mass loss) rather than the initial stellar mass as is the case in our definition of $m_{\mathrm{O}}$.

[TABLE]

as shown by the purple, dotted curves (identical in top and bottom panels). Both chemical evolution models assume instantaneous stellar mass loss and enrichment, with no time dependence (whereas our simulations consistently follow time-dependent mass loss and enrichment as the stellar population ages). The ratio $r/m_{\mathrm{O}}=26.7$ (using $r=0.4$, and $m_{\mathrm{O}}=0.015$, the fiducial values from Weinberg 2017) is thus the only free parameter. The models match the relative abundances at $z=0$ surprisingly well.

As can be seen by comparing the two panels of Figure 1, there is relatively little evolution in the correlation between deuterium and oxygen. At fixed oxygen metallicity, the deuterium abundance is slightly higher at $z=3$ than at $z=0$. This is because most of the oxygen is produced in core-collapse supernovae, which also dominate the stellar mass loss at early times. At late times, AGB stars are responsible for most of the mass loss, adding deuterium-free material, but not substantially enriching the gas with oxygen. This can be tested by dividing the cumulative amount of mass loss added to the gas in the simulations by the total amount of gas-phase oxygen at different redshifts. As mentioned before, the fiducial ratio used is $r/m_{\mathrm{O}}=26.7$, which is close to, though slightly higher than, the value we find at $z=0$, $r/m_{\mathrm{O}}=24.4$ (directly computed from and averaged over all our simulations). At $z=3$, however, the average ratio in our simulations is substantially different, $r/m_{\mathrm{O}}=19.4$. We therefore added extra model curves (cyan, dashed using Equation 4 and orange, dotted using Equation 5) to each panel of Figure 1, where we changed the value of $r/m_{\mathrm{O}}$ to match the average value in the simulations.

The level of agreement between these simple models, especially the closed box model from Tinsley (1980), and our cosmological simulation results is remarkable given the very different approaches. This lends credence to both methods and shows that the most important factor in this correlation is the ratio $r/m_{\mathrm{O}}$, which can be calculated from stellar population synthesis models. The complex processes involved in the formation of galaxies, such as galaxy mergers, time-variable star formation and galactic outflows, as well as the lack of mixing in these simulations are thus likely unimportant where these relative abundances are concerned. The improvement from the small variation in $r/m_{\mathrm{O}}$ with redshift supports our claim that the evolution is due to the extra (almost oxygen-free) mass loss from AGB stars at late times.

Accurate observations of (D/H) at $\mathrm{[O/H]}>-1$ in combination with an accurate determination of $\mathrm{(D/H)}_{\mathrm{prim}}$, either from observations at low metallicity or derived from CMB measurements, would be able to determine $r/m_{\mathrm{O}}$. The recycling fraction is governed by intermediate-mass stars as well as massive stars, whereas the oxygen yield depends only on the latter. Therefore, the relation between $\mathrm{(D/H)/(D/H)}_{\mathrm{prim}}$ and [O/H] can potentially be used to constrain stellar evolution models and/or the variation of the IMF at the high-mass end. Although all our simulations were run with the same stellar evolution model and IMF, numerical chemical evolution models have already demonstrated that the deuterium fraction depends on these choices (e.g. Tosi et al., 1998; Prantzos & Ishimaru, 2001; Romano et al., 2006).

In the top right panel, observational constraints independently derived by Linsky et al. (2006) (blue cross with error bars) and Prodanović et al. (2010) (red cross with error bars) for the local ISM have been included, slightly offset from $\mathrm{[O/H]}=0$ for clarity and using $\mathrm{(D/H)}_{\mathrm{prim}}=2.547\times 10^{-5}$, as recently obtained by Cooke et al. (2016). Our solar deuterium value is in excellent agreement with that of Linsky et al. (2006), but higher than that of Prodanović et al. (2010). The latter, however, is interpreted by the authors as a lower limit on the true value, in which case it is also in agreement with the result from our simulations (see Section 3.3 for further discussion on these results).

Additionally, it is of interest to compare our results to those from numerical chemical evolution models. We therefore reproduced one of the models, based on a Scalo 1986 IMF and Schaller et al. 1992 stellar lifetimes, from Romano et al. (2006), shown in the top right panel as the purple, dotted curve. Although these results agree qualitatively, there is a clear $2\sigma$ discrepancy at high metallicity. This is potentially caused by the different IMF and stellar evolution models used, by the different star formation histories (and thus different importance of AGB stars), or by the inclusion of galactic outflows. Given the excellent match between our simulations and the simple closed box model, we believe the former explanation is the most plausible.

Although the $1\sigma$ scatter in the relation between deuterium and oxygen abundances is small, it is nonzero and slightly larger at lower redshift. As mentioned before, almost all oxygen is produced by massive stars and released in core-collapse supernovae, whereas this is only the case for about half of the iron. The other half is synthesized in intermediate-mass stars and released in type Ia supernova explosions and winds from AGB stars (e.g. Wiersma et al., 2009). The iron abundance at fixed oxygen abundance therefore enables us to trace the relative importance of massive (younger) stars and intermediate-mass (older) stars and check whether variations in the contribution of stellar mass loss by AGB stars is responsible for the scatter seen in Figure 1.

Figure 2 shows the residual dependence of the normalized deuterium fraction on the iron metallicity at fixed oxygen metallicity at $z=0$ (top panel) and $z=3$ (bottom panel). $Z_{\mathrm{Fe}}/Z_{\mathrm{Fe,\astrosun}}$ is the iron mass fraction of the gas, normalized by the solar value. The coloured curves show the median relation between deuterium and iron in four bins, 0.1 dex wide, with (from top to bottom) $\mathrm{Log}_{10}Z_{\mathrm{O}}/Z_{\mathrm{O,\astrosun}}\approx-1$ (blue), $\approx-0.5$ (cyan), $\approx 0$ (black), and $\approx 0.5$ (red). Error bars show the 16th and 84th percentiles of the distribution and only bins containing at least 100 gas particles are included.

A larger iron abundance at fixed $Z_{\mathrm{O}}/Z_{\mathrm{O,\astrosun}}$ means that older, intermediate-mass stars have been relatively more important for enriching the ISM. The deuterium fraction is thus expected to decrease with increasing iron abundance. Figure 2 proves that this is indeed the case, although there is significant scatter in the residual relation between deuterium and iron. We therefore conclude that the (small) scatter in the relation between deuterium and oxygen is at least in part due to the varying importance of AGB stars and thus to the varying age of the stellar population responsible for enriching the gas.

Comparing the $z=0$ and $z=3$ results, it is clear that the residual dependence of the deuterium fraction on the iron abundance is stronger at $z=0$ than at $z=3$. This is likely because the variation in stellar population ages, and thus in the importance of AGB stars, is smaller at higher redshift, when the Universe was much younger. This is also consistent with the fact that the $1\sigma$ scatter in the correlation between the deuterium fraction and oxygen abundance is smaller at higher redshift. Measuring $\mathrm{[Fe/H]}$ besides $\mathrm{[O/H]}$ and (D/H) will provide even better constraints on the stellar IMF and stellar evolution models.

3.1.1 Milky Way-mass galaxies

To understand the chemical evolution of galaxies like the Milky Way, Figure 3 shows the deuterium evolution for three of our simulated galaxies with stellar masses close to that of the Milky Way at $z=0$, as indicated in the legend. The mean deuterium retention fraction is calculated for the gas within 20 proper kpc of the galaxies’ centres, with a temperature below $10^{4}$ K, which selects the warm ionized and cold neutral gas in the ISM. The black curves include the evolution of the hydrogen fraction, as it would be measured in observations. The orange curves show the recycled gas fraction, which is lower, but show the same trends with look-back time. The final $z=0$ values vary between the three galaxies, because they have different stellar masses. The mass dependence will be discussed in Section 3.3. Here, we are interested in the evolution of (D/H), that is, in the shape of the curves. Initially, the deuterium fraction is equal to its primordial value, after which it decreases. Two of the galaxies show an approximately linear decrease towards $z=0$ (solid and dashed curves), whereas for the galaxy with $M_{\mathrm{star}}=10^{10.4}$ M*☉* (dotted curve; ‘m12v’) the deuterium fraction levels off in the last $\approx 5$ Gyr. The former have therefore not reached an equilibrium between the inflow of deuterium-rich gas from the IGM, the addition of deuterium-free gas through stellar mass loss, and the outflow of deuterium-poor gas. The latter galaxy has potentially reached chemical equilibrium in its ISM.

Our three galaxies have different stellar masses and gas masses, on top of different star formation histories, and we lack the statistical power to control for this. Despite this limitation, we checked whether or not the low-redshift behaviour of the deuterium fraction is related to the star formation history. The galaxy which has reached deuterium equilibrium ($M_{\mathrm{star}}=10^{10.4}$ M*☉*; ‘m12v’) has already formed half of its stars by $z\approx 1.1$, whereas the other two reach half of their present-day stellar mass only at $z\approx 0.4$ and have thus experienced much more low-redshift star formation and thus more low-redshift stellar mass loss. Therefore, the reason for the different low-redshift behaviour may indeed lie in the different star formation histories of our simulated galaxies.

We also checked for a dependence of the deuterium evolution on the mass loading factor, i.e. the gas outflow rate from a galaxy divided by the galaxy’s SFR, as suggested by Weinberg (2017). The average mass loading factor at $z_{1}>z>z_{2}$ is calculated in the following way. We select all the gas particles within 20 proper kpc of the galaxies’ centres and a temperature below $10^{4}$ K. We then divide the total mass of those selected particles that have been turned into stars by $z=z_{2}$ by the total mass of the selected particles that are still gaseous, but located beyond 20 proper kpc of the galaxies’ centres at $z=z_{2}$. We take $z_{1}$ and $z_{2}$ to be approximately 1.5 Gyr apart, which is similar to the gas consumption time-scale. We find that the mass loading factor is relatively constant in the last 5 Gyr. From the most to least massive of our simulated Milky Way-mass galaxies, their average late-time mass loading factors over the last 5 Gyr are 0.2 (‘m12i’), 0.4 (‘m12v’), and 2.1 (‘m11.9a’). The values for ‘m12i’ and ‘m12v’ are consistent with the upper limits from Muratov et al. (2015), who argue that these low mass loading factors are not driven by galactic winds, but caused by random gas motions and/or close passages of satellite galaxies. We conclude that there is no clear correlation of the mass loading factor with the late-time deuterium evolution.

Knowing the evolution of (D/H) can thus potentially help us understand a galaxy’s star formation history. This could be achieved for the Milky Way with an accurate determination of the deuterium fraction in giant planets in the Solar System, such as Jupiter, in combination with present-day measurements in the local ISM (Lellouch et al., 2001). The deuterium fraction in the giant planets provides a fossil record of the deuterium fraction in the local ISM during the time the Solar System was formed, about 4.5 Gyr ago. Using $\mathrm{(D/H)}_{\mathrm{prim}}=2.547\times 10^{-5}$ from Cooke et al. (2016), the measurement by Lellouch et al. (2001) implies $\mathrm{(D/H)/(D/H)}_{\mathrm{prim}}=0.82^{+0.12}_{-0.15}$ in Jupiter, which is consistent with all three of our simulated Milky Way-mass galaxies and is not precise enough to distinguish between a declining or constant deuterium fraction. Future observations with higher accuracy would be well-suited for this purpose.

3.2 Deuterium fraction at high redshift

There has been a large observational effort to measure the deuterium fraction in metal-poor gas through absorption lines in spectra of background quasars. Lyman Limit Systems (LLSs; $10^{17.2}<N_{\mathrm{H\,\textsc{i}}}<10^{20.3}$ cm*-2*, where $N_{\mathrm{H\,\textsc{i}}}$ is the $H\,\textsc{i}$ column density) and Damped Lyman-$\alpha$ Systems (DLAs; $N_{\mathrm{H\,\textsc{i}}}>10^{20.3}$ cm*-2*) are optically thick to Lyman limit photons. To make a fair comparison between our simulations and these observations, we calculate column densities based only on the neutral gas. Because the gas comprising these strong absorbers is partially shielded from the ambient UV radiation, it is more neutral than if it were optically thin. This is taken into account in our simulations by using the fitting formula from Rahmati et al. (2013), which has been shown to capture the effect of self-shielding well.

Cooke et al. (2014) argue that the most precise measurements can be made in absorbers with $N_{\mathrm{H\,\textsc{i}}}>10^{19}$ cm*-2*. In order to compare to these systems, we also restrict ourselves to sightlines with column densities above this limit. Additionally, we discard the rare systems with $N_{\mathrm{H\,\textsc{i}}}>10^{21}$ cm*-2* in order to not be dominated by molecular gas. We note that neither this selection nor the self-shielding correction affects our results. We do not find a dependence of (D/H) on column density at fixed metallicity, so absorption line systems at any column density could be used. The vast majority of the selected high column density absorbers are located in the haloes around galaxies (van de Voort et al., 2012). We therefore use a simulated region of 300 by 300 proper kpc centred on the main galaxy in each of our zoom-in simulations. We grid this volume into 1 by 1 proper kpc pixels to calculate the column density of $\mathrm{H}\,\textsc{i}$, $\mathrm{D}\,\textsc{i}$, and $\mathrm{O}\,\textsc{i}$. We assume that the neutral fraction is the same for all three atoms, because their ionization potentials are very similar, as is also done in observations.

The black curve in Figure 4 shows the median fraction of deuterium in neutral gas, divided by its primordial value, in the selected LLSs and DLAs at $z=3$ as a function of their metallicity. The different grey scales show the $1\sigma$ (dark) and $3\sigma$ (light) scatter around the median. Observations of $(\mathrm{D}\,\textsc{i}/\mathrm{H}\,\textsc{i})$ compiled by Cooke et al. (2016) and their associated $1\sigma$ errors are shown as red error bars, where we assumed that $\mathrm{(D/H)}_{\mathrm{prim}}=2.547\times 10^{-5}$, the weighted mean of their measurements.

Our simulations confirm that at $\mathrm{[O/H]}\lesssim-2$ the deuterium abundance is very close to the primordial value (within 0.1 per cent), as seen before in Figure 1. These low-metallicity systems are therefore appropriate to use to determine $\mathrm{(D/H)}_{\mathrm{prim}}$. At $\mathrm{[O/H]}=-1$ the median deuterium abundance is still only 1 per cent below primordial, similar to the $1\sigma$ error in the weighted mean of the observational values from Cooke et al. (2016). The scatter in the relation is even smaller than in Figure 1, because we are including all (neutral) gas along a particular line-of-sight (rather than individual gas particles), decreasing the importance of small fluctuations. Observations of $N_{\mathrm{H\,\textsc{i}}}$, $N_{\mathrm{D\,\textsc{i}}}$, and $N_{\mathrm{O\,\textsc{i}}}$ are therefore well-suited to determine $\mathrm{(D/H)}_{\mathrm{prim}}$ and the relation between the deuterium and oxygen abundances.

Instead of assuming no variation as a function of metallicity for the 6 observed systems shown in Figure 4, we can test how well they match our simulations, which show a slight downward trend and minor additional scatter. We select those sightlines in our simulations that have the same metallicity as one of the observed absorbers, within $1\sigma$ errors. We then use least square fitting and calculate $\chi^{2}$ between our simulated sightlines and the observations as a function of $\mathrm{(D/H)}_{\mathrm{prim}}$, which sets the relative normalization. The minimum $\chi^{2}$ is reached for $\mathrm{(D/H)}_{\mathrm{prim}}=(2.549\pm 0.033)\times 10^{-5}$, where the errors are $1\sigma$ and calculated from the difference in $\chi^{2}$. This is consistent with theoretical models of Big Bang nucleosynthesis, based on cosmological parameters (Coc et al., 2015; Cyburt et al., 2016). Our best estimate is very similar to, though slightly higher than, the weighted mean calculated by Cooke et al. (2016), who assumed no metallicity dependence. For this low-metallicity sample, we do not gain much accuracy from comparing the data to our simulations.

However, given that the scatter in the simulations is much lower than the observational measurement error at all metallicities, more metal-rich absorption-line systems can be used to determine the primordial deuterium fraction. This would allow for the expansion of the observational sample, which would improve the accuracy of $\mathrm{(D/H)}_{\mathrm{prim}}$. Even absorbers with $\mathrm{[O/H]}\gtrsim-1$ can be used when taking into account the relation between the deuterium and oxygen abundance based on hydrodynamical simulations or on Equation 5 combined with a prescription for the change of $X_{\mathrm{gas}}$ with metallicity. For the latter, one should use a slowly evolving ratio of recycling fraction to oxygen yield, $r/m_{\mathrm{O}}$, increasing with time as the contribution of mass lost by intermediate-mass stars increases. Vice versa, observations of metal-rich absorbers can set constraints on the ratio of the recycling fraction and the oxygen yield, assuming that the primordial abundance of deuterium is known from either CMB measurements or from absorption-line observations at $\mathrm{[O/H]}\lesssim-2$. $r/m_{\mathrm{O}}$ depends on the relative number of intermediate- and high-mass stars and on their stellar yields and can thus potentially help constrain the high-mass end of the stellar IMF and/or stellar evolution models.

It is important to note that the depletion of deuterium onto dust and preferential incorporation into molecules could cause large scatter in (D/H) between quasar sightlines at fixed metallicity, which are not due to variations in the recycled gas fraction. This is probably seen in the local ISM at solar metallicity (e.g. Wood et al., 2004; Linsky et al., 2006; Prodanović et al., 2010) and briefly discussed in Section 3.3. Unfortunately we cannot address this issue with our current simulations. A relatively large sample of (D/H) measurements in absorption-line systems could quantify the scatter in (D/H) between sightlines at fixed metallicity. This will tell us whether the depletion of deuterium onto dust is important in the intergalactic medium at $\mathrm{[O/H]}>-2$, because our simulations have shown that the scatter due to variations in stellar mass loss at fixed metallicity is negligibly small. If dust depletion turns out to be dominant, these systems cannot be used for determining $\mathrm{(D/H)}_{\mathrm{prim}}$ or $r/m_{\mathrm{O}}$. However, there is some evidence that LLSs tend to reside in dust-poor environments (Fumagalli et al., 2016). Additionally, dust depletion seems to be less important in lower column density absorbers (e.g Linsky et al., 2006; Prodanović et al., 2010). It is therefore possible that these systems are well-suited for determining $\mathrm{(D/H)}_{\mathrm{prim}}$ or $r/m_{\mathrm{O}}$ even at $\mathrm{[O/H]}>-2$.

3.3 Deuterium fraction at low redshift

To compare with observations of (D/H) at $z=0$, we focus on the deuterium fraction in the ISM of our simulated Milky Way-like galaxies. The black curves in Figure 5 show how the ratio of the present-day abundance of deuterium to the primordial abundance varies with 3D distance from the galactic centre, $R_{\mathrm{GC}}$, for the same three galaxies as shown in Figure 3. For completeness, we show the recycled gas fraction as orange curves. $\mathrm{(D/D_{prim})}$ is similar to, though slightly lower than, $\mathrm{(D/H)/(D/H)}_{\mathrm{prim}}$, because of the decrease of the hydrogen fraction. One galaxy (dashed curve; ‘m11.9a’) has a central hole in its ISM, created by galactic winds, consistent with its relatively large average mass loading factor (see Section 3.1). It therefore has no deuterium measurement at $R_{\mathrm{GC}}<4$ kpc. This galaxy has the lowest deuterium abundance at large radii, because the gas that was originally in its centre has been moved to larger radii. The deuterium retention fraction is low in the centres of the other two galaxies, where the density of stars is high and most of the star formation takes place. The deuterium fraction for all three galaxies increases with galactocentric radius, as previously shown by chemical evolution models (e.g. Prantzos, 1996; Chiappini et al., 2002; Romano et al., 2006). The importance of stellar mass loss therefore increases towards the galaxy centre and recycled gas accounts for about half of the gas in the central kpc. The steepness of the deuterium abundance gradient could also reveal information on the assembly history of a galaxy (e.g. Prantzos, 1996). In our sample, the galaxy with the flattest deuterium profile has the highest outflow rate. It may therefore depend more strongly on (bursty) galactic outflows than on (smooth) gas accretion.

Large scatter exists between measurements of the deuterium abundance in the local ISM via absorption-line observations. This scatter could be explained by localized infall of pristine gas, with very little mixing. In this case, the average astration factor is relatively high (and mass lost from stars dominates the ISM). However, if this is the case, the oxygen abundance is also expected to decrease locally as it becomes diluted with the metal-free, infalling gas, resulting in large scatter. The fact that oxygen shows much smaller abundance variation than deuterium argues against such localized infall (Oliveira et al., 2005). Another, more likely, explanation for the large (D/H) sightline variations is that some of the deuterium is depleted onto dust. The probability of deuterium depletion onto dust grains and incorporation into molecules is high, since the zero-point energies of deuterium-metal bonds are lower than those of the corresponding hydrogen-metal bonds (Jura, 1982; Tielens, 1983). When the ISM is heated, dust grains and molecules can be destroyed, returning deuterium to the atomic gas phase. Metals, such as iron, silicon, and titanium, are also depleted onto dust grains and the correlation of their abundances with deuterium supports this theory (Prochaska et al., 2005; Linsky et al., 2006; Lallement et al., 2008). Based on the assumption that the observational scatter is caused by deuterium depletion onto dust, relatively high deuterium abundances, and low astration factors, are derived for the local Milky Way ISM by Linsky et al. (2006) and Prodanović et al. (2010).

The deuterium retention percentages in the solar neighbourhood, here defined as $7<R_{GC}<9$ kpc, lie between 85 and 92 per cent for our three simulations of star-forming galaxies with masses similar to that of the Milky Way444Although our calculations are done using 3D distance, the results are unchanged when we restrict ourselves to the ISM, because most of the gas mass lies within the star-forming disc.. This is consistent with Figure 1, where we found that 91 per cent of its primordial value is recovered for $\mathrm{[O/H]}=0$ at $z=0$. Using the value from Cooke et al. (2016) for the primordial deuterium fraction as in Section 3.2, the deuterium abundance derived by Linsky et al. (2006) implies that the local ISM still contains $91^{+9}_{-10}$ per cent of the primordial deuterium abundance. This is consistent with our simulations within $1\sigma$. Prodanović et al. (2010) use the same data compilation, but a different method, to derive a deuterium retention percentage of $79\pm 4$ per cent (again assuming $\mathrm{(D/H)}_{\mathrm{prim}}=2.547\times 10^{-5}$). This is consistent (within $2\sigma$) with our most massive Milky Way-like galaxy, with $M_{\mathrm{star}}=10^{10.8}$ M*☉*. However, Prodanović et al. (2010) stress that their measurement can also be interpreted as a lower limit in the event that all available sightlines are affected by dust depletion. In this case, our other galaxies are also consistent with their model. Our simulations exhibit low astration factors and therefore agree with the explanation that the large scatter in local ISM observations is due to dust depletion rather than due to poor mixing of freshly accreted gas.

No known galaxy besides the Milky Way has a measurement of the deuterium fraction in their ISM. Such observations would be interesting, because our simulations predict a strong dependence on stellar mass. Figure 6 shows the mean deuterium retention fraction (top panel) and recycled gas fraction (bottom panel) within 20 kpc of the centre of the galaxy for gas with a temperature below $10^{4}$ K as a function of stellar mass at $z=0$. Due to the depletion of hydrogen at higher metallicity, the differences between the deuterium retention fraction and recycled gas fraction increase with stellar mass, but the trends with mass remain the same. The black crosses show the mass-weighted mean, while the red diamonds show the (instantaneous) SFR-weighted mean. The latter is therefore a better indicator of how important stellar mass loss is for the fuelling of star formation, whereas the former is the value that would be measured, for example, in sightlines through the ISM. The galaxy shown in grey with the highest stellar mass has large uncertainties for its deuterium abundance, because it contains only 12 star-forming gas particles (but more than 1000 gas particles in total). Two other massive galaxies are not included, because they contained no star-forming gas and very little non-star-forming gas. Galaxies with $M_{\mathrm{star}}<10^{8}$ M*☉* are also excluded, because they contained no or very little star-forming gas. The other galaxies (shown in black, red, and blue) have 100 star-forming gas particles or more.

There is a definite trend with stellar mass, where more massive galaxies generally have lower deuterium fractions and thus a larger contribution of stellar mass loss. This is a clear prediction from our simulations. The exception is the lowest mass galaxy which has a low (D/H), as well as a high [O/H] (as expected from Figure 1). This reflects the variation in star formation, inflow, and outflow history. Our sample of galaxies is limited, but quantifying the scatter in the deuterium fraction between galaxies of similar stellar mass would help us to understand how this correlates with their formation history.

Due to their tight relation, the oxygen abundance could provide very similar information as the deuterium abundance. We expect this to work well at subsolar metallicities. At solar and supersolar metallicities, the dependence of (D/H) on [O/H] is very steep, which means that to determine the recycled gas fraction, very precise measurements of [O/H] are needed. Additionally, the scatter in the relation increases towards high metallicity at $z=0$, due to varying contributions of AGB stars, impeding the determination of the recycled gas fraction through [O/H].

The deuterium fraction is always lower when weighted by star formation rate (red diamonds) rather than by mass (black crosses), because this dense, star-forming gas is close to newly formed stars with high mass loss rates. This means that mass loss is more important for fuelling star formation than it is for replenishing the more diffuse ISM. A measurement of the deuterium fraction (such as those in the local ISM) generally provides information on the latter, but not the former.

In our simulations, the star formation in low-mass galaxies is fuelled predominantly by gas accretion, although there is a small contribution by stellar mass loss. For two of the galaxies with $M_{\mathrm{star}}\approx 10^{10.5}$ M*☉, mass loss fuels about half of the current star formation, while one other is still dominated by gas accretion. For the massive galaxy (or galaxies, if we include the one with large uncertainties) with $M_{\mathrm{star}}\approx 10^{11}$ M☉*, mass loss fuels the majority of star formation. As mentioned in Section 1, stellar mass loss has been suggested as fuel for the observed star formation in local, massive, early-type galaxies and in central cluster galaxies. Our simulations are consistent with this interpretation.

Notably, there are also detections of neutral deuterium in quasar sightlines through clouds outside the galactic disc of the Milky Way, in the lower galactic halo (Sembach et al., 2004; Savage et al., 2007). These clouds could provide the fuel to sustain the Milky Way’s steady star formation. However, the errors associated with these measurements are large and therefore not very constraining. More precise determinations of (D/H) could help constrain the nature of this gas. If the deuterium fraction of the gas around the Milky Way is high, it is likely pristine gas falling in from the IGM. If the deuterium fraction is similar to that of the Milky Way, it is likely part of a galactic fountain (Shapiro & Field, 1976). In intermediate cases, the gas could be a mix of both high and low (D/H) gas or the gas could be ejected or stripped from lower mass satellites (e.g. Anglés-Alcázar et al., 2017). These different theories can potentially be tested via estimates of the scatter in (D/H) in the halo, which would be large if the fountain gas is not fully mixed with the gas falling in from the IGM, but small if the gas was ejected or stripped from satellites.

In order to compare the deuterium fraction in the galaxy to the deuterium fraction of accreting gas, the blue triangles in Figure 6 show $\mathrm{(D/H)/(D/H)}_{\mathrm{prim}}$ at $z=0.1$ of the gas that is accreted at $0<z<0.1$, which means it is located at $R_{\mathrm{GC}}<20$ kpc and has $T<10^{4}$ K at $z=0$, but is located at larger radii and/or has higher temperatures at $z=0.1$. Our radial boundary is somewhat arbitrary, but our results are not very sensitive to this (except for the normalization, because the deuterium fraction is generally higher at larger radius). The accreting gas has a higher deuterium retention fraction than the ISM gas. This is consistent with a substantial part of the material accreting for the first time onto a galaxy. The fact that $\mathrm{(D/H)/(D/H)}_{\mathrm{prim}}<1$, however, shows that another part of the accreting gas has previously been ejected or stripped from a galaxy. Such a combination of different accretion channels has already been found in these and other simulations (e.g. Oppenheimer & Davé, 2008; Christensen et al., 2016; Anglés-Alcázar et al., 2017; van de Voort, 2017). The difference in (D/H) between accreting gas and the ISM is especially large for the more massive galaxies in our sample. As mentioned above, this balance between accretion from the IGM and reaccretion could potentially be shown observationally by measuring (D/H) in clouds that are accreting onto the Milky Way (e.g. Sembach et al., 2004).

4 Discussion and conclusions

We have quantified the evolution of the deuterium fraction and its dependence on the oxygen abundance, galactocentric radius, and stellar mass in cosmological zoom-in simulations with strong stellar feedback. Because deuterium is only synthesized in the early Universe, it provides an interesting constraint on cosmology and on galaxy evolution. The normalized deuterium fraction is a measure of the recycled gas fraction, i.e. the fractional contribution of stellar mass loss to the gas, because deuterium is completely destroyed in stars and therefore $f_{\mathrm{recycled}}=1-\mathrm{(D/D_{prim})}$. Observations, however, measure $\mathrm{(D/H)/(D/H)}_{\mathrm{prim}}$, which is higher than $\mathrm{(D/D_{prim})}$, especially at high metallicity (because $X_{\mathrm{gas}}/X_{\mathrm{prim}}\approx 1-3Z_{\mathrm{gas}}$). Our simulations self-consistently follow gas flows into and out of galaxies and the (metal-rich and deuterium-free) mass loss by supernovae and AGB stars. We have compared our predictions to available observations at low and high redshift and found them to be consistent. Our main conclusions can be summarized as follows:

• •

The deuterium fraction exhibits a tight correlation with the oxygen abundance, evolving slowly with redshift (Figure 1). This is captured well by simple chemical evolution models (Tinsley, 1980; Weinberg, 2017), which depend only on the ratio of the recycling fraction and the oxygen yield. We find a small increase in $r/m_{\mathrm{O}}$ with time, because of the increased importance of AGB stars. The variation in the importance of AGB stars at fixed oxygen abundance can be traced by the iron abundance and is responsible for some of the (small) scatter in the deuterium fraction (Figure 2).

• •

The three Milky Way-mass galaxies in our sample exhibit different evolution at low redshift (Figure 3). The galaxies that form many of their stars at late times have continually decreasing deuterium fractions in their ISM. The galaxy which forms most of its stars before $z=1$ shows a constant deuterium fraction in the last $\approx 5$ Gyr, indicating that it may have reached chemical equilibrium. The evolution of the deuterium fraction may therefore be directly connected to the galaxy’s star formation history.

• •

The deuterium fraction is very close to primordial at $\mathrm{[O/H]}\lesssim-2$ (within 0.1 per cent). These are the metallicities of LLSs and DLAs typically used to measure the primordial deuterium abundance. Because of the tight correlation with metallicity, deuterium measurements in more metal-rich systems can also be used to constrain the primordial deuterium fraction (Figure 4) if dust depletion is unimportant.

• •

We compared our simulations to the observational sample of Cooke et al. (2016) to determine a primordial deuterium fraction of $\mathrm{(D/H)}_{\mathrm{prim}}=(2.549\pm 0.033)\times 10^{-5}$, very close to, though slightly higher than, their original estimate, which assumed no dependence of the measured (D/H) on metallicity. Our result is also in agreement with cosmological parameters and Big Bang nucleosynthesis (Coc et al., 2015; Cyburt et al., 2016).

• •

The deuterium fraction increases with galactocentric radius. Our simulations are consistent with the available estimates from the local Milky Way ISM where the observed scatter between sightlines is assumed to be caused by depletion onto dust (Figure 5).

• •

The deuterium fraction decreases with increasing stellar mass, which means that the importance of stellar mass loss in our simulations increases with stellar mass (Figure 6). Mass loss is more important for fuelling star formation than for replenishing the general ISM (which has, on average, a lower gas density). Accreting gas has a higher deuterium fraction than the ISM of galaxies, but lower than primordial. This is consistent with previous findings that some gas accretes directly from the IGM, but some has been previously been ejected or stripped from a galaxy (e.g. Anglés-Alcázar et al., 2017).

Due to the tight correlation of (D/H) and [O/H], shown in Figure 1, measurements of the oxygen abundance could provide the same information as the deuterium abundance if this relation is accurately calibrated by observations at $\mathrm{[O/H]}\gtrsim-1$. The relation evolves slowly, because of the increased importance of stellar mass loss from AGB stars, which should be taken into account (see Section 3.1). If the importance of mass loss for fuelling the ISM increases for massive galaxies, as in our simulations, (D/H) will decrease with mass and [O/H] (and other metal abundances) will increase. There is observational evidence from gas-phase and stellar metallicities that this is indeed the case, although [O/H] may saturate at the highest stellar masses (e.g. Tremonti et al., 2004; Gallazzi et al., 2005; Mannucci et al., 2010; Peng et al., 2015). Using the median relation derived from our simulations (Figure 1) one can immediately estimate the contribution of stellar mass loss given the gas-phase oxygen metallicity. However, this only works in the situation where galactic winds remove mass loss from supernovae and AGB stars (approximately) equally, but may not if ejecta from young stars are removed from a galaxy (through supernovae or AGN, quenching star formation) after which its ISM is replenished by mass loss from old stars alone. This is likely the reason that the scatter in (D/H) increases at supersolar metallicity in our simulations. In observations, however, the scatter at high metallicity is probably dominated by the depletion of deuterium onto dust, an effect which is not included in our calculations.

Our cosmological, hydrodynamical simulations follow time-dependent chemical enrichment and the assembly of galaxies self-consistently, whereas simple, analytical chemical evolution models assume instantaneous recycling and specify a specific star formation history (Tinsley, 1980; Weinberg, 2017). We nevertheless compare our results to these models and find a remarkable agreement (especially with Tinsley 1980) when considering the relation between the deuterium and oxygen abundance, despite the very different methods used. Although other galaxy properties are not necessarily well reproduced in these simplified models (e.g. Binney & Merrifield, 1998), these do not play a dominant role when relative abundances are concerned. The bursty star formation and strong galactic outflows present in our simulations thus have no major impact on the correlation between (D/H) and [O/H]. Although Weinberg (2017) find that a relatively high mass loading factor is necessary to match the $z=0$ deuterium abundance in the local ISM, two of our simulated galaxies show low mass loading factors at late times (but high mass loading factors at early times, see Muratov et al. 2015).

Leitner & Kravtsov (2011) use zoom-in simulations with and without stellar mass loss to show that mass loss dominates the fuelling of star formation at late times for galaxies in haloes of similar mass to that of the Milky Way. However, their simulation did not include outflows from either star formation or AGN, which results in the galaxies being too massive at $z=0$. This means that too much baryonic mass is locked up in stars and therefore less gas is available for accretion at late times. Furthermore, the majority of the mass lost by stars is retained by the galaxy and not ejected by a galactic wind. Our simulations show that for Milky-Way mass galaxies with strong stellar feedback, cosmological inflow either dominates over or rivals stellar mass loss for the fuelling of star formation. However, recycled gas dominates the SFR at $M_{\mathrm{star}}\approx 10^{11}$ M*☉* in our simulations.

Segers et al. (2016) use large-volume simulations with stellar and AGN feedback that match the stellar–to–halo mass relation and find that the contribution of mass loss to the fuelling of star formation is largest for galaxies with $M_{\mathrm{star}}\approx 10^{10.5}$ M*☉*. The deuterium fraction in the ISM is therefore the lowest around this mass and increases for more massive galaxies. Our simulations (without AGN feedback) find a different behaviour at the high-mass end, where mass loss becomes increasingly important and the deuterium fraction decreases. Segers et al. (2016) also show results from simulations without AGN feedback, which do not match the stellar–to–halo mass relation, but agree qualitatively with our simulations. In the absence of AGN feedback, there are no strong galactic outflows, which means that most of the stellar mass loss is retained in the ISM. Different implementations of AGN feedback (ejective or preventive) could also result in different deuterium abundances. Future observations of (D/H) in massive galaxies therefore have the potential to discriminate between these different models.

Numerical chemical evolution models generally find somewhat lower deuterium fractions (higher astration factors) than we do in our cosmological simulations (e.g. Dearborn et al., 1996; Tosi et al., 1998; Chiappini et al., 2002; Romano et al., 2006). This mild difference could be due to the different assumptions made for the stellar IMF, metal yields, and stellar mass loss. A detailed quantitative comparison between the two different numerical approaches requires using the same IMF and stellar evolution models, which is left for future work. The star formation history of a galaxy determines the relative importance of AGB stars with respect to core-collapse supernovae, which may explain any remaining discrepancy. This could potentially be tested by comparing the oxygen and iron abundances (see Figure 2). Another possibility is that the inclusion of galactic outflows in our simulations reduces the importance of stellar mass loss, because a substantial fraction of the recycled gas is ejected from the ISM.

The large observed scatter in the deuterium fraction between different sightlines through the local ISM could be interpreted as evidence for an inhomogeneous ISM due to localized gas accretion, giving rise to a low average deuterium fraction (e.g. Hébrard & Moos, 2003). Our work instead favours the interpretation that the ISM is well-mixed and the observational scatter is caused by the depletion of deuterium onto dust, which leads to a local deuterium fraction not much lower than the primordial value (Linsky et al., 2006; Prodanović et al., 2010). It is also possible that both dust depletion and localized infall play a roll and the true value is intermediate between those cases (Steigman et al., 2007). More precise observations of the deuterium and metal abundances would help clarify which of these interpretations is correct.

In summary, we have quantified the deuterium fraction in a suite of zoom-in simulations and found it to be tightly correlated with the oxygen metallicity and consistent with current observational constraints. We conclude that the primordial deuterium fraction (and thus early cosmological expansion and Big Bang nucleosynthesis) can also be constrained by using observations at medium to high metallicity in combination with our simulations. Or, vice versa, if the primordial deuterium fraction is known, these measurements can inform us about the ratio of the recycling fraction to the oxygen yield and thus about the high-mass end of the stellar IMF and stellar evolution models. Our simulations predict that the deuterium fraction is lower at smaller galactocentric radii and for higher mass galaxies. This means that stellar mass loss could provide most of the fuel for star formation in massive early-type galaxies and in the centres of less massive, star-forming galaxies. Grid-based calculations or SPH simulations with explicit diffusion would be useful to determine whether or not small-scale mixing modifies the deuterium and oxygen abundances. Accurate observations of the deuterium fraction provide us with the possibility to understand the fuelling of star formation through stellar mass loss in galaxies in general and the Milky Way in particular.

Acknowledgements

We would like to thank the Simons Foundation and the organizers and participants of the Simons Symposium ‘Galactic Superwinds: Beyond Phenomenology’, in particular David Weinberg, for interesting discussions and inspiration for this work. We also thank Thomas Guillet and Joop Schaye for helpful discussions and Tim Davis for useful comments on an earlier version of the manuscript. We would like to thank the referees for valuable comments that helped clarify our results and put them into context. Support for FvdV was provided by the Klaus Tschira Foundation. EQ was supported by NASA ATP grant 12-APT12-0183, a Simons Investigator award from the Simons Foundation, and the David and Lucile Packard Foundation. CAFG was supported by NSF through grants AST-1412836 and AST- 1517491 and by NASA through grant NNX15AB22G. DK was supported by the NSF through grant AST-1412153 and by the Cottrell Scholar Award from the Research Corporation for Science Advancement. Support for PFH was provided by an Alfred P. Sloan Research Fellowship, NASA ATP Grant NNX14AH35G, and NSF Collaborative Research Grant #1411920 and CAREER grant #1455342. Numerical calculations were run on the Caltech compute cluster “Zwicky” (NSF MRI award #PHY-0960291), through allocation TG-AST120025, TG-AST130039 and TG-AST150045 granted by the Extreme Science and Engineering Discovery Environment (XSEDE) supported by the NSF, and through NASA High-End Computing (HEC) allocation SMD-14-5189, SMD-14-5492, SMD-15-5950, and SMD-16-7592 provided by the NASA Advanced Supercomputing (NAS) Division at Ames Research Center.