THE HISTORY OF R-PROCESS ENRICHMENT IN THE MILKY WAY

We use the high resolution, zoom-in cosmological simulation of a Milky Way Galaxy analog "Eris" to track the production and transportation of r-process elements. A detailed description of the Eris simulation is provided by Guedes et al. (2011). Here we briefly outline the aspects relevant to this study. The simulation was performed with the parallel TreeSPH code Gasoline (Wadsley et al. 2004) in a WMAP-3 cosmology. The run includes a uniform UV background (Haardt & Madau 1996), Compton cooling, atomic cooling, and metallicity-dependent radiative cooling at T $\lt \;{10}^{4}$ K. Star formation is modelled by stochastically forming "star particles" out of gas that is sufficiently cold ($T\lt 3\times {10}^{4}$ K) and reaches a threshold density of ${n}_{\mathrm{SF}}=5$ atoms cm⁻³. The local star formation rate (SFR) follows $d{\rho }_{*}/{dt}=0.1{\rho }_{\mathrm{gas}}/{t}_{\mathrm{dyn}}\propto {\rho }_{\mathrm{gas}}^{1.5}$, where ${\rho }_{*}$ and ${\rho }_{\mathrm{gas}}$ are the stellar and gas densities, respectively, and ${t}_{\mathrm{dyn}}$ is the local dynamical time. Each star particle has initial mass ${m}_{*}=6000\ {M}_{\odot }$ and represents a simple stellar population that follows a Kroupa et al. (1993) initial mass function (IMF), and inherits the metallicity of its parent gas particle. Star particles inject energy, mass, and metals back into the ISM through SNe Ia, SNe II, and stellar winds (Stinson et al. 2006). Eris' high resolution enables the development of an inhomogeneous ISM which allows realistic clustered star formation and strong cumulative feedback from coeval supernova explosions. Large-scale galactic winds are launched as a consequence of stellar feedback, which transports a substantial quantity of metals into the circumgalactic medium and enriches the subsequent gas accretion (Shen et al. 2013). At z = 0, Eris forms an extended, rotationally supported stellar disk with a small bulge-to-disk ratio. The structural properties, the mass budget in various components, and the scaling relations in Eris are simultaneously consistent with observations of the Galaxy (Guedes et al. 2011).

The simulation follows Raiteri et al. (1996) to model metal enrichment from SNe II and SNe Ia. Metals are distributed to gas within the smoothed particle hydrodynamics (SPH) smoothing kernel (which consists of 32 neighboring particles). For SNe II, metals are released as the main-sequence progenitors die, and iron and oxygen are produced according to the following fits to the Woosley & Weaver (1995) yields:

$$\begin{eqnarray}&&{M}_{\mathrm{Fe}}=2.802\times {10}^{-4}{\left(\displaystyle \frac{{m}_{*}}{{M}_{\odot }}\right)}^{1.864}{M}_{\odot },\end{eqnarray} \tag{ 1 }$$

and

$$\begin{eqnarray}&&{M}_{{\rm{O}}}=4.586\times {10}^{-4}{\left(\displaystyle \frac{{m}_{*}}{{M}_{\odot }}\right)}^{2.721}{M}_{\odot }.\end{eqnarray} \tag{ 2 }$$

For SNe Ia, each explosion produces 0.63 ${M}_{\odot }$ of iron and 0.13 ${M}_{\odot }$ of oxygen (Thielemann et al. 1986). Stellar wind feedback is based on Kennicutt et al. (1994), and the returned mass fraction was determined following Weidemann (1987). The returned gas inherits the metallicity of the star particle. We adopt the Asplund et al. (2009) solar abundance scale for elements other than O and Fe, which are not tracked in the simulation.

We note that there are several limitations to our modeling. First, because the smallest gas resolution element is an ensemble of gas particles within the smoothing kernel rather than individual particles, forcing the newly formed star particle to inherit metallicity only from its parent gas particle may amplify sub-resolution metallicity variances between gas particles. Second, the simulation used a traditional SPH formalism where metals advect with the fluid perfectly, without mixing due to microscopic motions (such as turbulence). Both caveats may cause an artificially inhomogeneous chemical distribution (Wiersma et al. 2009; Shen et al. 2010). While improved runs within the Eris suite include a model for turbulent diffusion (Shen et al. 2013), for this study we have used a simulation without mixing so that our r-process injection method is consistent with the production and distribution of O and Fe in the simulation. In Section 3, we explore a simple diffusion model to illustrate the effect of mixing. Instead of using the metallicity of the parent gas particle, we average the metallicity over 128 neighboring gas particles to determine the metallicity of the newly formed star particle.

In addition, despite being one of the highest-resolution galaxy formation simulations that evolves to z = 0, our model still lacks the spatial and temporal resolution to correctly follow the formation and enrichment of the first generation of Population III stars as well as early Population II stars. It is a common practice in the literature to introduce a metallicity "floor" at high redshift to account for the unresolved earliest population (e.g., Krumholz & Gnedin 2011; Kuhlen et al. 2013; Hopkins et al. 2014). The metallicity floor is typically around ${10}^{-4}$–${10}^{-3}\;{Z}_{\odot }$, motivated by detailed simulations of Population III star formation and primordial metal enrichment (e.g., Wise et al. 2012). Although a metallicity floor was not included in Eris during the run, in Section 3 we investigate the effect of a modest metallicity floor (${10}^{-4}\;{Z}_{\odot }$) implemented in post-processing. This is achieved by assigning a minimum abundance of [Fe/H] = −4.0 to every gas particle, and an α-enhancement of [O/Fe] + 0.4, corresponding to the IMF-weighted alpha enhancement estimated from models of zero-metallicity stellar nucleosynthesis (Woosley & Weaver 1995; Heger & Woosley 2010; Limongi & Chieffi 2012).

The key ingredient in our analysis is the realistic star formation history (SFH) of a galaxy simulated in a cosmological context, which at redshift z = 0 is a close analog of the Milky Way (Guedes et al. 2011). In this section we describe our post-processing implementation for NS mergers. The important elements of our model include: (1) the delay-time distribution (DTD) of mergers; (2) the merger rate and the yield of r-process elements; and (3) the spatial distribution of injection sites and their sphere of influence.

2.2.1. Delay-time Distribution, Merger Rates, and R-process Yields

The merger DTD, P(t), is well modeled by a power law (Piran 1992; Kalogera et al. 2001), although there is significant uncertainty on the value of the exponent (e.g., Behroozi et al. 2014, and references therein). Herein, we adopt $P(t)\propto {t}^{-n}$ for $t\gt {t}_{\mathrm{cut}}$ and zero probability otherwise, where ${t}_{\mathrm{cut}}$ is the initial time delay after a burst of star formation before the first merger occurs. Herein we consider two power-law indices, n = 2 and n = 1, and assume conservative values of ${t}_{\mathrm{cut}}=100$ Myr and ${t}_{\mathrm{cut}}=200$ Myr (Fryer et al. 1999; Belczynski et al. 2006).

Our fiducial model assumes that each merger event produces a mass ${M}_{\mathrm{rp}}=0.05\;{M}_{\odot }$ of r-process elements (e.g., Just et al. 2015). We also consider a model where ${M}_{\mathrm{rp}}=0.01\;{M}_{\odot }$ is produced, which reflects the mass of dynamically ejected material (Lattimer & Schramm 1974; Rosswog et al. 1999; Metzger et al. 2010; Roberts et al. 2011; Bauswein et al. 2013; Grossman et al. 2014; Ramirez-Ruiz et al. 2015). For each merger event, we track the production of europium, which we assume is produced in solar relative proportions such that ${M}_{\mathrm{Eu}}/{M}_{\mathrm{rp}}=9.3\times {10}^{-4}$ (Sneden et al. 2008).

The NS merger rate is then determined by convolving the SFH extracted from the simulation with the DTD:

$$\begin{eqnarray}&&{\mathcal{R}}(t)=A{\int }_{{t}_{\mathrm{cut}}}^{{t}_{{\rm{H}}}}{\dot{M}}_{*}(t-\tau )\;P(\tau )\;d\tau \end{eqnarray} \tag{ 3 }$$

where A is a constant that is fixed by the total number of merger events, ${\dot{M}}_{*}$ is the SFR, and t_H is the Hubble time. To calculate A, we impose that the abundance of Eu/O in the final simulation output corresponds to the solar value (i.e., [Eu/O] = 0.0). There are ∼260 million massive stars formed in Eris that end their life as a SN II, and the IMF-weighted O yield per SN II event is $\simeq 1\ {M}_{\odot }$ in our model. Assuming a solar number ratio, log(Eu/O)${}_{\odot }$ = −8.14, the total mass of Eu needed to obtain the solar [Eu/O] abundance is $\sim 18\;{M}_{\odot }$. Thus, the total mass of r-process elements produced during the chemical evolution of Eris needs to be $18,800\;{M}_{\odot }$. If we now assume that each compact binary merger contributes ${M}_{\mathrm{rp}}=0.05\;{M}_{\odot }(0.01\;{M}_{\odot })$ of r-process, then we require $\approx 3.76\times {10}^{5}(1.88\times {10}^{6})$ mergers in 13.8 Gyr to explain the observed solar Eu/O ratio. The constant A is therefore fixed by requiring that the integral of the merger rate (Equation (3)) over the lifetime of Eris is equal to $3.76\times {10}^{5}\;(1.88\times {10}^{6})$. The resulting NS merger history is shown in Figure 1 as the solid red curve, for ${M}_{\mathrm{rp}}=0.05\;{M}_{\odot }$, which is in good agreement with the expected rates calculated by Abadie et al. (2010). For reference, we also present the Eris SN II rate as a function of time as the blue curve in Figure 1.

Zoom In Zoom Out Reset image size

2.2.2. Injection History

Compact binary mergers are expected to predominantly occur within ∼10–100 kpc of a Milky Way-like host galaxy (Bloom et al. 1999; Belczynski et al. 2006; Kelley et al. 2010). Therefore, the spatial distribution of NS mergers broadly follows the stellar distribution of the host galaxy. We have thus adopted a post-processing implementation to include NS mergers in Eris. Our approach is justified since the momentum imparted to the surrounding gas by a merger is much less than that of a SN explosion, despite releasing an energy that is similar in magnitude to a SN. The gas dynamics is therefore largely unchanged.

In the top left panel of Figure 2, we present a face-on illustration of the projected surface mass density of Eris' stars at redshift z = 2. By construction, this panel represents the distribution of merger injection sites. Similarly, in the top right panel of this figure, we present the SFR surface density, which traces SN II. In the bottom panels of Figure 2, we show the surface mass density of gas on two different spatial scales, which are enriched by these events and later form a new generation of stars.

Zoom In Zoom Out Reset image size

Using the above formalism, we calculate the number of NS mergers that occur between adjacent timesteps and randomly select a corresponding number of star particles from a uniform distribution to act as the merger injection sites. The surrounding gas particles are then enriched with a total Eu mass ${M}_{\mathrm{Eu}}^{\mathrm{tot}}$ = $4.65\times {10}^{-5}\;{M}_{\odot }$ (corresponding to ${M}_{\mathrm{rp}}$ = 0.05 ${M}_{\odot }$), which is distributed over the 32 neighboring gas particles according to the smoothing kernel, as outlined in Wadsley et al. (2004). We note that the oxygen and iron enrichment follow an identical scheme. The evolution of each gas particle is tracked for subsequent timesteps. When a stellar particle is born in the non-diffusion case, it inherits the Eu/Fe ratio from the parent gas particle. In our simple mixing model, the stellar particle inherits the average abundance of Eu, O, and Fe from 128 neighboring gas particles as described in Section 2.1. Hereafter, when comparing with observations, we only consider star particles present in the z = 0 snapshot.

In our investigation, we have found that compact binary mergers are able to successfully produce a large dispersion in the [Eu/Fe] ratios at low metallicity, which is diluted at higher metallicity due to chemical mixing. In general, the large r-process yield and the rarity of NS mergers can produce a highly varied enrichment of the r-process. We have also explored extensions to our fiducial model including: (1) a DTD with $P(t)\propto {t}^{-2}$ (top left panel of Figure 5); (2) increasing the frequency of mergers by a factor of five, while reducing the yield per event by the same amount (top right panel of Figure 5); (3) doubling the initial time delay for mergers after a burst of star formation, to ${t}_{\mathrm{cut}}=200\;\mathrm{Myr}$ (bottom left panel of Figure 5); and (4) adopting a less efficient mixing algorithm (bottom right panel of Figure 5). A general result of decreasing the power-law index of the DTD is that more events will occur at earlier times (i.e., shortly after the 100 Myr delay that we have imposed). The elevated number of events at early times acts to modestly reduce the dispersion in [Eu/Fe] at low metallicity. Similarly, by decreasing the Eu yield per event, the number of events is globally increased in order to match the observed Eu/O ratio in the Sun. The increased number of events also acts to slightly reduce the dispersion in [Eu/Fe] at all metallicities. When we change the merger delay time to be ${t}_{\mathrm{cut}}=200\;\mathrm{Myr}$, the evolution of [Eu/Fe] is qualitatively very similar to our fiducial model for [Fe/H] $\gtrsim -2.0$. However, at early times and hence lower metallicities ([Fe/H] ≲ –2.0), there are subtle differences between this model extension and the fiducial model, which are largely because there are fewer lower metallicity stars being enriched with Eu. Finally, if we reduce the efficiency of mixing relative to the fiducial model, by averaging the metallicity of the 32 nearest gas particles when a star particle is born, our model provides a poorer description of the [Eu/Fe] chemical evolution, particularly for the low [Eu/Fe] stars.⁵

Zoom In Zoom Out Reset image size

We therefore conclude that modest variations of the model parameters for the injection history of the r-process do not significantly alter our results—Eris displays a general correspondence with the chemical evolution of [Eu/Fe] observed in the Milky Way. The choice of mixing does not significantly affect the [Eu/Fe] distribution of the strongly Eu-enhanced stars, where [Eu/Fe] ≳ 0.0, and both models provide an acceptable fit to the data in this regime. The choice of mixing does, however, affect the distribution of the lowest [Eu/Fe] stars. This regime is also the most poorly constrained by observations, since there are many stars where the Eu abundance has not been measured (e.g., Barklem et al. 2005; François et al. 2007; Cohen et al. 2008; Roederer et al. 2010, 2014). The level of mixing is therefore best determined by using the abundance distribution of other elements, such as the [α/Fe] ratio as implemented herein. We encourage future studies that consider a wider range of chemical elements to obtain a better handle on the mixing process.

The chemical inhomogeneity of the interstellar medium is a key aspect of chemical evolution that can be better addressed in numerical simulations. We now explore the effect of chemical inhomogeneities and metal mixing in more detail. In particular, we consider a scenario that is similar to a 1D chemical evolution model where the metal content of all star-forming gas (i.e., gas that turns into stars at the next simulation timestep) is completely mixed before forming a new generation of stars. This assumption has been used in many analytic and semianalytic calculations of chemical evolution. Such models provide a greater degree of physical intuition, but are unable to capture all of the relevant physical processes that are important for GCE; for example, these models are not typically set in a cosmological setting. The results of our 1D experiment are shown in Figure 6, where the dark and light contours display the 68% and 95% contours enclosing the chemical evolution of our fiducial model, and the single blue line illustrates the result for our 1D experiment. At late times, the centroids of the two distributions agree very well, although there is clearly a very large dispersion in all of the element abundances that are considered in this work. However, at early times, and hence lower metallicities, our completely mixed chemical evolution model provides a very poor resemblance to the fiducial model. Moreover, the 1D model is unable to keep track of the dispersion at early times (see also Argast et al. 2004), which is particularly important for Eu.

Zoom In Zoom Out Reset image size

Figure 7 further illustrates the shortcomings of 1D chemical evolution models by plotting the evolution of O and Eu relative to Fe. It is evident that in the 1D models, stars with metallicities [Fe/H] ≲ −3.0 cannot be formed after ≳100 Myr and that most of the star-forming material in the Milky Way has reached a level of [Fe/H] ∼ −3.0 in the z ∼ 15–20 range. The ubiquity of Eu in stars with [Fe/H] ≲ −3.0 in our cosmological simulations indicates that metal diffusion is not widespread from the sites of NS mergers until after a few Gyr. This heavy metal flow, mixed and sheared by gas motions, remains highly anisotropic and is responsible for generating the larger-scale Eu abundance variations pervading low [Fe/H] stars. The red curve in Figure 7 contrasts the model Mod1NS' from Matteucci et al. (2014) with our cosmological chemical evolution. Their Mod1NS' model shares the closest similarity to our NS merger prescription and as such provides the clearest model comparison. It is clear that 1D chemical evolution models may misrepresent the true chemical evolution of a galaxy, particularly at early times when metal mixing is not yet effective. Future studies of detailed chemical evolution in a cosmological setting (e.g., Kobayashi & Nakasato 2011; Maio et al. 2015) are thus encouraged to study the properties of metal-poor gas and stars.

Zoom In Zoom Out Reset image size

After the submission of this paper we became aware of a recent preprint (van de Voort et al. 2015), in which the authors simultaneously reached similar conclusions to those presented herein regarding the role of NS mergers in Galactic r-process enrichment using the cosmological zoom-in FIRE simulation (Hopkins et al. 2014). From the parameter study by van de Voort et al. (2015), they also note that the abundance ratio and scatter of the r-process/Fe ratio at low metallicity may be sensitive to uncertainties in modeling turbulent mixing of the ISM in their simulation. In Figure 8, we compare the chemical evolution of Eu/Fe for their fiducial model relative to ours. Both simulations appear to track the solar abundance of Eu/Fe for [Fe/H] ≳ −2.0, and exhibit a similar level of scatter. Although both simulations exhibit a large dispersion of [Eu/Fe] values for [Fe/H] ≲ −2.0, the centroid of the two simulations shows an opposite evolution, which likely reflects the different model choices of Eu injection sites, numerical resolution, SFH, and feedback prescriptions between the two simulations.

Zoom In Zoom Out Reset image size

Like all chemical evolution models, our conclusions may be sensitive to the assumptions of our model, such as the implementation of metal diffusion as well as the other numerical choices outlined in Section 2.1. It is nevertheless encouraging that our chemodynamical simulation with a simple, first-order prescription of subgrid metal mixing provides an acceptable agreement to the observed chemical evolution of [α/Fe]. Future cosmological simulations that follow the chemical evolution of elements produced on different timescales and from a larger variety of production sites may help to elucidate some of the finer details of chemical mixing.

Our results demonstrate the importance of tracing both the dynamics and chemical evolution self-consistently in studying the origins and the evolution of r-process elements. In this regard, chemodynamic simulations are advantageous in modeling the chemical inhomogeneity in the ISM, and are highly complementary to computationally inexpensive analytic calculations that provide a somewhat higher level of physical intuition. Future analytic and semianalytic calculations that include a prescription of the effects of chemical inhomogeneity and mixing may provide a very useful insight into galactic chemical evolution, especially for trace elements, such as the r-process.

While there need not be a single r-process production site, the current data favor a source that is either uncommon or whose yield only sporadically produced the r-process, as argued by many other authors in the past. Our study has found that NS mergers are a strong candidate as a dominant production site, even at low metallicity. Future cosmological hydrodynamic simulations with a more realistic mixing prescription, which follows a galaxy with a similar chemical evolution and SFH to that experienced by the Milky Way, are now required to investigate this problem in further detail.