The Origin of r-process Elements in the Milky Way

Fits to the UV, optical, and IR data from GW170817 provide definitive constraints on the masses and general composition of ejecta components. Table 1 summarizes the relevant findings in recent literature. All studies agree that at least two spatially separated ejecta components were present: low-opacity radioactive material to power the early optical emission (Evans et al. 2017; Nicholl et al. 2017), and high-opacity lanthanide-polluted outflow to account for late-time near-IR emission (Arcavi et al. 2017; Tanvir et al. 2017; Troja et al. 2017). Note that only the latter contains Eu, while the former component is lanthanide-free and expected to have a relatively high Y_e.

Studies with bolometric light curve reconstruction (Cowperthwaite et al. 2017; Kasliwal et al. 2017b; Rosswog et al. 2017; Smartt et al. 2017) provide simple and robust lower limits on the total mass of the ejecta (m_total > 0.03–0.05 M_⊙) and upper limits on the neutron richness, measured in terms of electron fraction: Y_e ≲ 0.3 (Rosswog et al. 2017). The latter constraint follows from the shape of the nuclear heating profile, as reflected in the bolometric light curve.

For better agreement with observations, some studies argued that an additional third ejecta component is required, such as a wide-angle mildly relativistic cocoon (Kasliwal et al. 2017b), secular disk winds (Perego et al. 2017), or a "purple kilonova" outflow (Cowperthwaite et al. 2017). All these additional outflows have relatively high electron fractions and do not produce Eu. In a different approach, an axisymmetric two-component model with toroidal dynamical ejecta and spherical wind was applied (Wollaeger et al. 2017). Geometric effects in this model allow photon reprocessing and thus higher luminosities for the same masses, correspondingly producing lower-mass estimates (see Table 1). In particular, based on this model, Tanvir et al. (2017) and Troja et al. (2017) required a moderate amount of high-Y_e wind ejecta (∼0.015 M_⊙), with a relatively low amount of low-Y_e dynamical ejecta (∼0.002–0.01 M_⊙). The statistical MCMC analysis of a large set of gray-opacity models in Cowperthwaite et al. (2017) and Perego et al. (2017), on the other hand, argued for a total ejecta mass of 0.04 M_⊙, with 1% of this mass in lanthanides. The mass of the low-Y_e component can also be estimated independently from fits of synthetic spectra at late times (Chornock et al. 2017; Kasliwal et al. 2017b; Tanaka et al. 2017).

In this paper, we test the hypothesis that GW170817 is a typical, representative event that produces a regular r-process signature, consistent with the robust abundance pattern observed in metal-poor halo stars (Sneden et al. 2008) and the ultra-faint dwarf galaxy Reticulum II (Ji et al. 2016; Roederer et al. 2016). This case is studied in Section 4. The other possibility of GW170817 being an unusual, non-representative case is discussed in Section 5.4, where theoretical r-process yields and their uncertainties (see Section 2.3) are applied. From these abundance patterns, we extract Eu yields in order to calibrate the NS–NS merger rates required by GCE studies, as they typically use Eu to trace the r-process enrichment (see Section 5.5).

Based on the estimates presented in Table 1, we adopt a conservatively broad range of 0.002–0.01 M_⊙ for dynamical ejecta, and 0.01–0.03 M_⊙ for the high-Y_e component, which adds up to 0.01–0.04 M_⊙ for the total ejecta mass range. Table 2 summarizes the resulting yields of r-process constituents. The second column displays the observed composition of r-process residuals in the solar system (from Arnould et al. 2007). The third column contains a hypothetical normalized composition expected if the low-Y_e component robustly produces the observed pattern (Korobkin et al. 2012; Bauswein et al. 2013; Mendoza-Temis et al. 2015), starting from the second peak on (A > 110). The fourth and fifth columns contain theoretical compositions for low-Y_e and high-Y_e components, computed from first principles with a nucleosynthesis network (see Section 2.3 for details). In the last two columns, we convolve the ejecta mass ranges inferred from GW170817 with the observed or theoretical mass fractions of different constituents. The penultimate column (${M}_{\mathrm{ejected}}^{\mathrm{obs}}$) lists masses computed assuming the observed solar r-process residuals, while the last column (${M}_{\mathrm{ejected}}^{{\rm{w}}/\mathrm{nuc}.}$) lists the ranges computed theoretically.

Table 2.  Yields of the r-process Constituents

Abundance	Xobs [10−8]	${X}_{\mathrm{low}\,{Y}_{e}}^{\mathrm{obs}}$	${X}_{\mathrm{low}\,{Y}_{e}}^{{\rm{w}}/\mathrm{nuc}.}$	${X}_{\mathrm{high}\,{Y}_{e}}^{{\rm{w}}/\mathrm{nuc}.}$	${M}_{\mathrm{ejected}}^{\mathrm{obs}}$ [M⊙]	Mejectedw/nuc. [M⊙]
Total r-process (A > 79)	${35.0}_{-0.3}^{+0.4}$	0.99a	0.98–1.0a	0.56–0.70	0.01–0.04	0.0075–0.03
Main r-process (A > 130)	${7.08}_{-0.03}^{+0.09}$	${0.284}_{-0.003}^{+0.003}$	0.70–0.99	0.002–0.06	(6–30) × 10−4	0.0014–0.012
1st peak (30 < Z < 38)	${7.92}_{-0.04}^{+0.09}$	0.0	0.0–2 × 10−4	0.7 − 0.9	0.004–0.012	0.007–0.03
2nd peak (48 < Z < 59)	${3.50}_{-0.05}^{+0.07}$	${0.141}_{-0.002}^{+0.002}$	0.2–0.8	0.007–0.1	0.004–0.013	0.0005–0.011
3rd peak (74 < Z < 83)	${1.62}_{-0.03}^{+0.20}$	${0.065}_{-0.001}^{+0.006}$	0.13–0.5	0.0	(1.3–7) × 10−4	(0.3–5) × 10−3
Trans-lead (Z > 82)	${0.03}_{-0.02}^{+0.06}$	${1.2}_{-0.6}^{+1.6}\times {10}^{-3}$	0.006–0.154	0.0	(2–30) × 10−6	(0.12–15) × 10−4
Iron Peak (21 < Z < 30)	0.0	0.0	0.0	0.0002–0.01	0.0	(5–14) × 10−6
Europium (Z = 63)	${0.036}_{-0.001}^{+0.005}$	${1.45}_{-0.003}^{+0.14}\times {10}^{-3}$	0.0002–0.02	0.0	(3–15) × 10−6	(0.4–22) × 10−6
56Ni	0.0	0.0	0.0	0.0	0.0	0.0

Note.  Second column: mass fractions X^obs in the solar r-process residuals (Arnould et al. 2007). Third column: mass fractions in the hypothetical low-Y_e dynamical ejecta when adopting the r-process residual abundance pattern for A > 110. Fourth and fifth columns: mass fractions with theoretical nucleosynthetic yields with low-Y_e and high-Y_e ejecta, respectively. The two remaining columns show the ejected masses inferred for GW170817: ${M}_{\mathrm{ejected}}^{\mathrm{obs}}$—using the observed r-process residuals, ${M}_{\mathrm{ejected}}^{{\rm{w}}/\mathrm{nuc}.}$—using theoretical yields. A and Z represent the mass and atomic numbers, respectively.

^aThe remaining ≈1% of mass is theorized to be in ⁴He formed by α-decays.

Download table as:  ASCIITypeset image

Note that the high-Y_e "wind" ejecta observed in GW170817 does not contain lanthanides and thus is not expected to produce third-peak r-process elements. In principle, the ratio of third to second r-process peaks from GW170817 may not match the solar abundance pattern. The uncertainties in wind versus dynamical ejecta masses for current light curve fits to GW170817 shown in Table 1 allow for a wide range of ratios of the r-process peaks (see Section 5.5 for discussion). However, this does not clash with the robust abundance pattern observed in the r-process-rich Galactic halo stars and in Reticulum II (Ji et al. 2016), because the robustness only holds for Z ≥ 56 (see e.g., Figure 11 in Sneden et al. 2008). It is currently unclear whether the robustness can be extended to the second r-process peak, since atomic lines in the region 50 ≤ Z ≤ 55 are in the UV range, which is notoriously hard to measure (e.g., Roederer & Lawler 2012; Roederer et al. 2012). In fact, some fission models of the robust r-process nucleosynthesis in the neutron-rich ejecta strongly underproduce the second peak (Goriely et al. 2013).

The consistency between the r-process pattern observed in metal-poor halo stars (Sneden et al. 2008), the solar r-process residuals (Arnould et al. 2007), and the ultra-faint dwarf galaxy Reticulum II (Ji et al. 2016; Roederer et al. 2016) places tight constraints on the site of the main r-process. To agree with observations, a feasible candidate needs to be able to reproduce the observed pattern robustly, with little sensitivity to the variations in system parameters (Korobkin et al. 2012). However, the current nuclear theory of heavy element nucleosynthesis produces uncertainties in the predicted pattern exceeding observational constraints by at least one order of magnitude (Mumpower et al. 2016b). In this section, we explore the approach of extracting the r-process yields of GW170817 using nucleosynthesis calculations from first principles.

The nuclear physics uncertainties mainly stem from the fact that the r-process path meanders through the uncharted territory of heavy, extremely neutron-rich nuclei close to the neutron drip line, for which no experimental data are available. Variations in the unknown nuclear masses (Mumpower et al. 2016b), fission fragment distribution (e.g., Goriely et al. 2013), neutron capture rates (Mumpower et al. 2012), β-decay rates (Mumpower et al. 2014; Eichler et al. 2015), and the specifics of the fission mechanism itself can all significantly impact nucleosynthetic yields (Mumpower et al. 2016b). For this reason, we calculate abundance yields with a variety of nuclear mass models and fission fragment distributions.

Mass model uncertainties are of utmost concern for nucleosynthesis calculations, as they influence all major features of the r-process abundance pattern. In particular, the strengths of the shell closures predicted by each model directly determines the shape of the primary r-process peaks, while more subtle trends for highly deformed nuclei influence the rare-earth peak. We select 10 popular mass models (DZ, FRDM1995, FRDM2012, HFB17, HFB21, HFB24, WS3, KTUY, SLY4, and UNEDF0) that are based on a wide range of physical underpinnings.

Some models, such as FRDM and WS3, are based on the macroscopic-microscopic approach that separates the macroscopic shape degrees of freedom from the microscopic description of the nucleus. Other models, such as SLY4 or UNEDF0, provide a fully microscopic description of the nucleus based on Skyrme interactions and density functional theory. The predictions of nuclear masses between these models tends to diverge as one approaches the neutron drip line resulting in a range of abundance patterns in calculations of nucleosynthesis.

Since nuclear binding (mass) is critical input for other nuclear properties relevant to the r-process, we self-consistently update nuclear reactions (Mumpower et al. 2017) and decay modes (Mumpower et al. 2016a) as in Mumpower et al. (2015) using the Los Alamos suite of statistical Hauser-Feshbach codes (Kawano 2016). We also include β-delayed and neutron-induced fission rates calculated within this same self-consistent framework. To accommodate the flexibility in nuclear input data required for these self-consistent calculations, we use the nucleosynthesis code PRISM (Portable Routines for Integrated nucleoSynthesis Modeling, Mumpower et al. 2017). This code easily allows for all nuclear data to be user-defined and permits a straightforward hierarchy so that theoretical values are overwritten with experimental ones (e.g., from NUBASE Audi et al. 2017), when available.

We first consider the astrophysical conditions of a wind scenario with Y_e = 0.27 previously found to produce the best fit to a "blue" kilonova light curve (Evans et al. 2017; Wollaeger et al. 2017). The bands in Figure 1 demonstrate the range of abundance predictions when different mass models are considered. Even with an order of magnitude variation in the dynamical timescale (different band colors), this environment is not capable of generating a full r-process. Only the second peak r-process elements are reached and Eu is not produced.

Zoom In Zoom Out Reset image size

Since the electromagnetic counterpart of GW170817 was found to be consistent with lanthanide production, we next turn to neutron-rich dynamical ejecta and consider two different astrophysical conditions. The first one is a "cold" merger outflow condition, such as what can be found in the tidal tail ejecta, with an entropy of 10 k_B/baryon and Y_e = 0.05 (Just et al. 2015). The second one is an astrophysical trajectory for "slow" ejecta (Mendoza-Temis et al. 2015) with a neutron-to-seed ratio of R_n/s ∼ 10³, which takes into account the reheating due to nuclear reactions with a temperature rise from ∼0.2 GK to ∼1.5 GK in ∼60 ms for the case considered here.

The "cold" trajectory dynamics are such that the temperature and density drops make most of the r-process proceed under an interplay between neutron capture and β-decay, since photodissociation falls out of equilibrium early. However, a scenario with reheating extends the time in which $(n,\gamma )\leftrightarrow (\gamma ,n)$ equilibrium persists allowing for more late-time neutron capture that shifts and narrows the third peak relative to the the results in "cold" conditions (see the different colored bands in Figure 2). The bands in Figure 2 demonstrate that the sensitivity of reaction rates to nuclear mass can alone produce an order of magnitude variation in the abundance for the considered dynamical ejecta trajectories.

Zoom In Zoom Out Reset image size

Additional complications arise in these very neutron-rich environments since the r-process path proceeds to the region of the nuclear chart where fission reactions dominate, and different fission prescriptions lead to distinct abundance patterns. The left and right panels of Figure 2 compare abundance predictions using the fission fragment distribution of Kodama & Takahashi (1975) versus assuming a simple, symmetric split for the fissioning nucleus. The fission prescription sets the distribution of the second r-process peak and the left edge of the rare-earth peak. The predictions from all mass models follow the common trend determined by the fission fragment distribution. On average, variations due to the fission prescription can lead to an overproduction (left panel) or underproduction (right panel) of Eu (Z = 63).

For GW170817, the constraints on the second and third peak r-process yields are set by the amount of high-opacity lanthanide elements needed to explain the late-time "red" kilonova emission. Uncertainties in the nuclear cross-sections can produce the same total lanthanide ejecta, but vary the production of individual components wildly. Keeping the amount of total lanthanides equal (their mass is constrained by the observations), we can study the additional uncertainty in the trans-lead, Eu, and r-process peak elements. The final column of Table 2 shows the yield range, including the nuclear physics uncertainties outlined in Table 3, which are based on the range of abundance predictions given by the nuclear mass models and fission fragment distributions outlined in this section. In this range of models, we find that the total Eu abundance can increase or decrease by a factor of 2 with the same total lanthanide abundance.

Population synthesis models predict NS–NS merger rates for stellar populations (e.g., Fryer et al. 1999; Voss & Tauris 2003; Dominik et al. 2012; Mennekens & Vanbeveren 2014; Belczynski et al. 2016a; Chruslinska et al. 2018). Those models can be confronted with the observed merger rate estimated from several known NS–NS systems in the Milky Way (${21}_{-14}^{+28}$ Myr⁻¹; Kim et al. 2015). For comparison with other observational constraints such as short-duration gamma-ray bursts (Berger 2014) and gravitational wave measurements (Abbott et al. 2016), a calculation of cosmological NS–NS merger rate densities is required. This involves tracing the formation of NS–NS progenitor systems according to the cosmic star formation history (CSFH, Madau & Dickinson 2014) and following their evolution until they merge using metallicity-dependent delay-time distributions (DTDs; see Belczynski et al. 2016b).

The merger rate densities based on previous calculations (Belczynski et al. 2016a) are too low compared to the latest LIGO/Virgo estimates at low redshift⁹ (${1540}_{-1220}^{+3200}$ Gpc⁻³ yr⁻¹; Abbott et al. 2017b). Those models have been revisited by Chruslinska et al. (2018). For many realizations of the input physics (e.g., varied assumptions about the natal kicks, angular momentum loss, mass transfer), the classical evolution of isolated binaries typically leads to low merger rate densities at low redshifts (≲50 Gpc⁻³ yr⁻¹). However, several models with specific common envelope physics, low angular momentum loss during Roche-lobe overflow, electron-capture SNe allowed in a wide range of initial stellar masses (with no natal kick applied), and reduced natal kicks for NS progenitors with heavily stripped envelopes, can produce local NS–NS merger rate densities as high as ∼500–600 Gpc⁻³ yr⁻¹.

Uncertainties associated with the CSFH, the stellar initial mass function, the binary fraction, and the evolution of metallicity through cosmic time can further shift the predicted merger rate densities by a factor of ∼2. The highest merger rate density predicted with the calculations of Chruslinska et al. (2018) is shown as the upper limit of the green shaded area in Figure 3. This limit (∼10³ Gpc⁻³ yr⁻¹ at redshift z = 0) represents the most optimistic model increased by a factor of 2 to show the currently attainable maximum NS–NS merger rate density with population synthesis methods. The lower limit is the same as in C17a (but see Chruslinska et al. 2018).

Zoom In Zoom Out Reset image size

Since the DTD of NS–NS mergers typically follows a power law in the form of t⁻¹ (see also Dominik et al. 2012), close NS–NS binaries merge within a few Myr after their formation. The evolution of the merger rate density should therefore follow the CSFH, which peaks at z ∼ 2. However, because NS–NS systems are most efficiently formed at high metallicities at z < 2, the peak of the merger rate density is shifted to z ∼ 1.5.

We use a similar approach to calculate the merger rate densities required by GCE. However, instead of using the NS–NS merger rates predicted by population synthesis for individual stellar populations, we use the ones adopted in GCE simulations. Those rates are calibrated to reproduce the [Eu/Fe]¹⁰ abundances observed in the Milky Way, assuming NS–NS mergers are the only source of r-process elements. We chose this particular elemental ratio for probing the r-process production because it represents the common observational target used in the seven GCE studies compiled and normalized in C17a.

The merger rate density required to reproduce the current [Eu/Fe] abundances in our Galaxy depends on the DTD of NS–NS mergers, on the chemical evolution code, and on the amount of Eu and Fe ejected by NS–NS mergers and supernovae, respectively (see Sections 5.3 and 7 in C17a for more details). The range of solutions for a DTD in the form of t⁻¹ is shown as the dark and light blue shaded areas in Figure 3 (see Section 4 for details).

These rates, however, do not account for other sources of uncertainties such as the rate of Fe injection by Type Ia supernovae (e.g., Matteucci et al. 2009; Côté et al. 2016; Rybizki et al. 2017), the fraction of NS–NS binaries merging outside the star-forming region (e.g., Fryer et al. 1999; Behroozi et al. 2014; Safarzadeh & Côté 2017, but see Beniamini et al. 2016), and the CSFH, which could all increase the width of the blue uncertainty bands.

Instead of using GCE simulations and reproducing stellar abundances, analytical calculations can also be used to estimate the role of NS–NS mergers on the Galactic r-process enrichment (e.g., Metzger et al. 2010; Bauswein et al. 2013; Rosswog et al. 2013; Hotokezaka et al. 2015). Those calculations are based on the total mass of r-process elements M_r,tot currently present the Milky Way, which can be obtained by multiplying the total mass of baryons in our Galaxy by the mass fraction of the r-process in the solar system (Qian 2000). This quantity can be divided by the lifetime of the Galaxy (∼10 Gyr) to calculate the average mass injection rate of r-process events. Then, by adopting different ejected masses for NS–NS mergers, one can infer the merger rates needed to recover M_r,tot.

By combining the properties of GW170817 with analytical calculations similar to the ones described above, several studies have demonstrated that NS–NS mergers can synthesize enough r-process material to be the dominant site in the Milky Way (Abbott et al. 2017a; Chornock et al. 2017; Cowperthwaite et al. 2017; Gompertz et al. 2017; Kasen et al. 2017; Rosswog et al. 2017; Tanaka et al. 2017; Wang et al. 2017; Hotokezaka et al. 2018). GCE simulations and analytical calculations are thus converging toward a similar message. However, uncertainties in the mass ejected by NS–NS mergers are still significantly affecting the predictive power of both approaches.

In this work, we used the conservative mass range of 0.002–0.01 M_⊙ for the dynamical ejecta of GW170817, which represents the low-Y_e component that synthesizes Eu (see Section 2.1). However, as seen in Table 1, some studies derived larger values that can reach up to ∼0.04 M_⊙. Using this larger value yields 4 times more Eu compared to our upper limit of 1.5 × 10⁻⁵ M_⊙. In that case, the NS–NS merger rate density required by GCE drops to ∼100 Gpc⁻³ yr⁻¹, which is about 3 times lower than the lower limit established by LIGO/Virgo. Therefore, an upper limit of 0.01 M_⊙ for the dynamical ejecta seems to be more consistent with the current observed merger rate (see the dark blue shaded area in Figure 3).

However, as discussed in Section 5.4, we cannot exclude the possibility that GW170817 is an unusual event, in which case the merger rate density and average mass ejected per NS–NS merger event would be subject to changes in the upcoming years.

Although population synthesis models overlap with GCE and LIGO/Virgo (see the green shaded area in Figure 3), their predicted merger rate densities tend to be lower. As a matter of fact, when addressing NS–NS merger rates in the context of old spheroidal galaxies with extinct star formation, similar to NGC 4993, which hosted GW170817, population synthesis models predict rates that are significantly lower than the ones established by LIGO/Virgo (Belczynski et al. 2017). This discrepancy could be solved by improving the physics behind the formation of NS–NS mergers in theoretical calculations, or by favoring mergers with longer delay times. The latter option would, however, potentially be in tension with the chemical evolution trend of Eu in the Milky Way, which overall seems to behave like an alpha element with short production timescales (e.g., Battistini & Bensby 2016; Spina et al. 2018, see also discussion in C17a and Hotokezaka et al. 2018).

So far, only one NS–NS has been detected by LIGO/Virgo, which means that GW170817 could be an unusual and non-representative event. The merger rate density could actually be lower, meaning that GW170817 has been detected earlier than statistically expected. If this is the case, the derived NS–NS merger rate will decrease as the LIGO/Virgo's observing time gets longer. From a GCE perspective, given the uncertainties, NS–NS mergers could still be the main r-process site even if the merger rate density was reduced, as long as it does not drop below ∼100–200 Gpc⁻³ yr⁻¹. On the other hand, if GW170817 has been detected later than statistically expected, the actual merger rate density could be higher.

GW170817 could also be unusual in terms of its total mass ejected. If NS–NS mergers eject on average more or less mass than GW170817, the range of merger rate densities required by GCE would be modified. Indeed, depending on the masses and mass ratio of the two neutron stars in the merger, the dynamical and wind ejecta masses can vary by a factor of 2–4 (Korobkin et al. 2012). The uncertain inferred neutron star mass ratios range from 0.4 to 1.0 (Abbott et al. 2017b). If the mass ratio is closer to 0.4, the ejected mass may be higher than representative values.

To connect GCE with LIGO/Virgo's detection, we used Eu as the r-process tracer to calibrate the required NS–NS merger rates. We recall that in GCE simulations, only one r-process abundance pattern is typically applied to all NS–NS mergers. To define whether these events are frequent enough to explain the r-process in the Milky Way, the adopted abundance pattern needs to represent the average yields synthesized by NS–NS mergers. In that regard, it is justified to use the solar r-process residuals.

However, as mentioned in Section 2.2, it is not guaranteed that the third to second r-process peak ratio ejected by GW170817 follows the solar distribution. In the case where the second peak is overestimated relative to the third peak, using Eu from Table 2 with GCE to quantify the role of NS–NS mergers would be irrelevant, as matching Eu would lead to an overestimation of the second peak. On the other hand, if the second peak in the GW170817 ejecta is underestimated, matching Eu with GCE would be reliable, but additional r-process sites would be needed to generate the missing lighter r-process elements (e.g., neutrino-driven winds in CC SNe, Arcones & Thielemann 2013).

With only one confirmed merger event, it is difficult to establish what is the typical r-process ejecta of NS–NS mergers. Since mergers can produce a variety of r-process abundance patterns, one would need to constrain the probability distributions of the ejecta before testing the contribution of NS–NS mergers using GCE simulations. Theoretical calculations could be used to define such distributions, but variations caused by nuclear physics uncertainties first need to be reduced. From an observational point of view, there are variations in the r-process abundance patterns of metal-poor stars when comparing with the light and heavy parts of the solar r-process composition (e.g., Sneden et al. 2008; Roederer et al. 2010; Hansen et al. 2014). Whether variations in NS–NS merger yields could participate in this observational feature needs further investigation.

As for GCE, we note that analytical calculations (see Section 5.1) also currently rely on the solar r-process residuals to test whether NS–NS mergers are the dominant r-process site.