Possibility to Identify the Contributions from Collapsars, Supernovae, and Neutron Star Mergers from the Evolution of the r-process Mass Abundance Distribution

The origin by rapid neutron capture (r-process) of nearly half of the heavy atomic nuclides from iron to uranium remains an open question (Kajino et al. 2019). The neutrino-driven wind (NDW) of core-collapse supernovae (CCSNe) is now believed to at best produce only light r-process elements (Wanajo 2013). The magnetohydrodynamically driven jet (MHDJ) from rapidly rotating, strongly magnetized CCSNe is an alternative site (Winteler et al. 2012). Also, there is a growing consensus (Frebel 2018; Hotokezaka et al. 2018) that neutron star mergers (NSMs) could be the source of r-process elements in the galaxy. This is due in part to the discovery of gravitational waves from the NSM GW170817 and its associated kilonova and GRB 170817A (Abbott et al. 2017a, 2017b).

In addition to the above sources, a single massive star collapsing to a black hole (collapsar) may also be a viable site for the main r-process abundances (Siegel et al. 2019). The definition of the term "collapsar" is ambiguous, but in this paper we take it to refer to a star that collapses gravitationally and becomes a black hole without leaving a neutron star. The type of associated explosion depends on the physical properties of the progenitor. If there is enough ejected material, it can be observed as a superluminous supernova (SLSN), while if almost all of the surrounding material falls into the black hole, it would be a dim transient called a failed supernova.

Indeed, it has been argued (Côté et al. 2018) that an extra production site for the r-process may be required. In particular, the data seem to require a source that was active in the early galaxy but fades away at higher metallicity. Otherwise, it is not possible to account for the decrease of r-process elements with iron after the onset of Type Ia supernovae (at [Fe/H] ∼ −1).

The universal elemental abundance pattern observed in metal-poor halo stars and the solar system (Sneden et al. 2008) suggests that only a single site contributed to the r-process elements (Mathews et al. 1992; Argast et al. 2004; Ishimaru et al. 2005; Côté et al. 2018; Hotokezaka et al. 2018). However, the isotopic abundance patterns from various sources can be quite different, even while maintaining the universal behavior of the elemental Z-distributions (Shibagaki et al. 2016; Suzuki et al. 2018; Siegel et al. 2019). This motivates the study described herein of the relative contributions of multiple r-process sites (CCSNe, NSMs, and collapsars) and their cosmic evolution.

There is a fundamental difference between NSMs, CCSNe, and collapsars. There is an unavoidable time delay from formation of the progenitors until the ejection of r-process material. CCSNe (producing an NDW and/or MHDJ) and collapsars result when a single massive star completes its evolution within a few Myr. They can enrich r-process elements in the interstellar medium (ISM) of the Galactic disk from the earliest times.

On the other hand, NSMs involve the remnants of previously exploded massive stars. Also, the observed orbital properties of binary pulsars imply a coalescence timescale ranging from a few hundred Myr to longer than the Hubble time (Swiggum et al. 2015). Even in the limit that the shortest timescale for mergers is only the stellar evolution timescale (∼10⁶ yr), it is inevitable that neutron star binaries are formed with a distribution of separation distances and merger timescales (e.g., Mathews et al. 1992; Argast et al. 2004). This leads to a longer delay in the arrival of r-process material in the solar neighborhood.

A new feature emphasized in this study is that the isotopic mass distribution can be used to test the source of r-process material at low metallicity. Indeed, the next generation of large telescopes equipped with a high-resolution spectrograph (Skidmore et al. 2015) may be able to provide this kind of information.

To illustrate our point we have utilized a widely employed galactic chemical evolution (GCE) model (Timmes et al. 1995) adapted to calculate r-process contributions from various sources in the solar neighborhood. We have considered some reasonable possible variations in the parameter space of r-process models, nuclear inputs, and astrophysical parameters to examine the general features of the evolution of the r-process mass distributions. We then impose the constraint that the universality of the elemental distribution should be reproduced.

We find as others have done that CCSNe and collapsars dominate the r-process abundances in the early galaxy, while NSMs can only arrive later when the metallicity has already been enriched to −1 < [Fe/H]. We show that this is independent of the minimum timescale for binary neutron-star coalescence and is a consequence of the unavoidable distribution of binary separation distances at formation. Moreover, we note that there could be unique features in the mass distribution function at low-metallicity. These features could be used to confirm the relative contribution of various sources to the r-process at low metallicity.

In the next Section 2, we formulate our GCE model including the NSM r-process as well as CCSN and collapsar r-process. We briefly review the population synthesis of the binary systems and define the coalescence time delay for NSM. In Section 3, we discuss the time evolution of both the isotopic mass (A) distribution and the elemental (Z) abundance patterns of the r-process nuclei that we calculated in our GCE model. In Section 4, we discuss the relative contributions to the r-process nuclei from CCSNe, collapsars, and NSMs. We also show the differences among the abundance evolution of typical r- and s-process elements, and in Section 5, we summarize the conclusions of the present work.

The observed cosmic star formation rate (SFR; Madau & Dickinson 2014) as well as galactic chemodynamical evolution (GCDE) models using smoothed particle hydrodynamics (SPH) of spiral galaxy evolution (Kobayashi 2004) and dwarf spheroidal galaxy evolution (Hirai et al. 2015) all indicate that the SFR at first rises and then diminishes with time. The GCE model (Timmes et al. 1995) adopted in this study produces an SFR similar to that deduced from observations and GCDE simulations.

The adopted GCE model also reproduces well the chemical evolution of light elements from hydrogen to zinc (Timmes et al. 1995). In this study we have extended this model to include the r-process contributions from NDWs, MHDJs, collapsars, and NSMs.

In this model, the gas evolution involves a cycle of star formation, stellar evolution, and nucleosynthesis; ejection of material into the interstellar medium (ISM); mixing of ejecta with the ISM; and formation of the next generations of stars. We adopt an exponentially declining galactic inflow rate with a timescale of 4 Gyr consistent with the hierarchical clustering paradigm. Although, the merger of dwarf galaxies into the galactic halo can bring some r-process-enriched stars, the bulk of the inflowing gas consists of pristine r-process-depleted material from the circumgalactic medium (Péroux & Howk 2020).

The surface density σ_i of isotope i in the ISM then obeys

$$\begin{eqnarray}\begin{array}{rcl}{\dot{\sigma }}_{i}(t) & = & \ \displaystyle \sum _{\mu }\ {\epsilon }_{\mu }\ {\displaystyle \int }_{{m}_{l}}^{{m}_{h}}\ {E}_{i,\mu }\ B(t-\tau (m))\ \phi (m)\ {dm}\\ & & +\,{\epsilon }_{\mathrm{NSM}}\ {\displaystyle \iiint }_{{M}_{l}}^{{M}_{h}}{dq}\ {da}\ {P}_{a}(a)\ {P}_{q}(q)\\ & & \times {E}_{i,\mathrm{NSM}}\ B(t-\tau ({m}_{2})-{\tau }_{g}(a))\ \phi (M)\ {dM}\ \\ & & -\,B(t)\ \displaystyle \frac{{\sigma }_{i}}{\ {\sigma }_{\mathrm{gas}}\ },\end{array}\end{eqnarray} \tag{ 1 }$$

where the first and second terms on the right-hand side describe the enrichment of newly produced nuclei by explosive nucleosynthesis in different astrophysical sites, i.e., μ =NDWs, MHDJs, collapsars, and NSMs, respectively. The third term accounts for the loss from the ISM due to star formation, where σ_gas is the total gas surface density. E_i,μ is the yield of isotope i from each astrophysical site μ. The quantity B(t) is the star formation rate, and ϕ(m) is the initial mass function that we adopt from Kroupa (2001).

For the present illustration we adopt the abundance distribution of r-process nuclei in the NDW model from the 1.8M_⊙ proto-neutron star yields of Wanajo (2013). We adopt yields from the NSM models of Suzuki et al. (2018) and Shibagaki et al. (2016), and the MHDJ model yields are taken from Nishimura et al. (2012). For the collapsar model, we adopt the yields of Nakamura et al. (2015), Famiano et al. (2020), and Siegel et al. (2019).

In Equation (1), the enrichment of r-process nuclei from supernovae and collapsars in the ISM is delayed by the period τ(m) from star formation to the death of the progenitor star of mass m. We adopt lifetimes from Schaller et al. (1992) for massive stars and from Woosley & Weaver (1995) for stars with m ≤ 10 M_⊙. In Woosley & Weaver (1995) it was estimated that the metallicity effect on the progenitor lifetime is only about ∼5% and is ignored.

The quantities _μ and _NSM are the fractions of stars resulting in each event in a given mass range m_l ∼ m_h or M_l ∼ M_h for a single or binary system. In this study, we take these as adjustable efficiency parameters, which are fit to the observed event rates of each site at the present time so as to account for the later solar r-process abundances. We here define the abundance of solar r-process nuclei as the sum of nuclides in the mass range 90 ≤ A ≤ 209.

In the present study we do not consider the slight metallicity dependence of various contributions. For the NSM contribution in Equation (1) we include the long gravitational-wave coalescence time τ_g(a) in addition to τ(m), where a is an initial binary separation. We adopt a paradigm whereby main-sequence binaries of total mass M = m₁ + m₂ are formed within a gas cloud. The heavier star with mass m₁ explodes first to form a neutron star or black hole, followed by the second CCSN of the lighter progenitor.

Binary population synthesis studies indicate that the minimum coalescence time is ∼100 Myr. This is in reasonable agreement with the lower limit estimated from observed binary pulsars (Swiggum et al. 2015). In the quadrupole formula for gravitational radiation (Peters & Mathews 1963), the coalescence time τ_g calculated from an initial binary separation a scales as a⁴. We ignore the dependence on eccentricity as tidal interactions should circularize the orbits.

The event rate of NSMs in Equation (1) includes the uniform probability P_q (q) for a given mass ratio of the secondary to total mass, q ≡ m₂/M, and the probability P_a (a) ∝ 1/a of an initial separation a. This is consistent with the observationally inferred coalescence time, ${P}_{{\tau }_{g}}\propto {\tau }_{g}^{-1}$ (Beniamini & Piran 2019; Simonetti et al. 2019). The lower limit of the initial separation a then constrains the minimum coalescence time.

The observed lower limit on the coalescence timescale of binary neutron stars is ∼50 Myr (Swiggum et al. 2015). However, there are several effects that might shorten the coalescence time. If the binary system is perturbed by a third stellar object or if there exists a common envelope, these would accelerate orbital energy-loss rate. We therefore treated the minimum value of a as a free parameter, and allowed shorter minimum coalesce times, i.e., $\min ({\tau }_{{\rm{g}}})=1$, 10, and 100 Myr in our GCE calculations. Figures shown in later sections are for the case of this value equal to 10 Myr. This does not change any of the general features of the GCE. Indeed, the results do not change even if we increase to the observed minimum value of 50 Myr.

We take efficiency parameters _μ in Equation (1) as free parameters adjusted so as to best fit the event rate of each site at the present time. The event rate of supernovae (SNe) in the galaxy is well studied. Recently, the frequency of NSMs has also been estimated by the detection of gravitational waves (GWs) and/or short gamma-ray bursts (SGRBs) (Hotokezaka et al. 2018).

Models for MHDJs and collapsars have large uncertainties. Both of them are considered as a source of SLSNe and long gamma-ray bursts (LGRBs). In this paper we consider two combinations of the event rate of each astrophysical site for r-process nucleosynthesis: One set of models (Model 1) assumes that collapsars are failed supernovae. For this we adopt a modified version (Famiano et al. 2020) of the collapsar r-process model of Nakamura et al. (2015). In this case MHDJ SNe are the source of SLSNe. In the second model (Model 2) collapsars are taken to be SLSNe. For this we adopt the r-process yields of Siegel et al. (2019).

Differences in the time delay from star formation to nucleosynthesis lead to variations in the isotopic r-process abundance patterns with metallicity. In general, the delay is shortest for collapsars and longest for NSMs. Therefore, galaxies in the early universe have a large contribution from collapsars, while the contribution from the NSMs is largest at the present time.

Figures 1 and 2 show examples of the abundance mass distribution as a function of metallicity and for different fission fragment distributions (FFDs) for both Model 1 (Figure 1) and Model 2 (Figure 2). The lines on these figures correspond to the NSM contribution (red line), collapsar contribution (blue line), and the CCSNe contribution (green line). The fit of the total abundances to the solar system r-process abundances is shown by a black line.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Among the features to be noted on these figures, the NSM and collapsar abundances have the possibility to avoid the underproduction near A = 140–145 that appears in many previous r-process calculations (Kajino et al. 2019). This is because the r-process occurs in an extremely neutron-rich environment. Thus, it undergoes fission recycling so many times that the fission fragments can enhance the abundances at A ≈ 130–150.

The tidal component of dynamical ejecta in NSMs has a very low electron fraction (Y_e ∼ 0.1). This causes the r-process path to run along extremely neutron-rich isotopes. The reaction flow quickly reaches heavy fissile nuclei in the region A = 250–290. Therefore, the FFD strongly affects the final abundance distribution in the mass range of A = 100–180 and even heavier isotopes (Shibagaki et al. 2016). When a symmetric FFD (Suzuki et al. 2018) is used in the NSM r-process, the second peak around A ∼ 130 is reproduced reasonably well as shown in the left panels of Figures 1 and 2, while the abundance distribution is smoothed out in the model with an asymmetric FFD (Shibagaki et al. 2016) in the right panels. This potentially solves the long-standing problem in r-process nucleosynthesis of underproduction of nuclides with A =140–145.

Nevertheless, there are several common features independent of the choice of FFDs; the third peak is slightly shifted toward the heavier mass region since the neutron number density remains high after the freeze-out in tidal ejecta from NSMs. The abundance patterns of the collapsars have not only these features but also another feature caused by their high entropy environment; an odd–even pattern manifests in the lanthanide abundance hill. This is a typical feature resulting from very rapid neutron captures and the sudden freeze-out. We note, however, that the calculation of Siegel et al. (2019) adopts a different ejecta model and shows a smoother pattern (Figure 2).

Although NSMs produce r-process contributions from not only the tidal ejecta, but also from accretion disk outflows as neutrino-driven winds (Kajino et al. 2019), we focus here on the dynamical tidal ejecta as in Suzuki et al. (2018) with 30 outflow trajectories based upon the hydrodynamic simulations of Korobkin et al. (2012), Piran et al. (2013), and Rosswog et al. (2013). These trajectories are based on SPH simulations in Newtonian gravity, where the neutrino transport is taken into account in a neutrino leakage scheme (Rosswog & Liebendoerfer 2003).

There is growing evidence in numerical simulations that both wind and jet outflow components contribute to the r-process nucleosynthesis (Martin et al. 2015). These components strongly affect the lighter r-process elements including the first peak, but they do not significantly produce the heavy elements above the second and third peaks and the hill of rare-earth elements (Thielemann et al. 2017). Actually, they could change the abundance pattern for the nuclei having A < 120 in our Figures 1–3. However, their impact would be negligibly small in our model since the NDW from CCSNe is the dominant contributor to this lower-mass region (A < 120).

Zoom In Zoom Out Reset image size

The properties of neutron-rich unstable nuclei and the physical conditions of the trajectory dynamics are thus intricately connected in any nucleosynthesis simulations. However, if some characteristic features, such as those discussed here, were measured in the isotopic abundances of r-process-enhanced metal-deficient stars, this could confirm that the heavy nuclei in those stars originate from a single or a few very similar r-process events.

3.1. The Variation of Isotopic Abundances

3.2. Universality in Elemental Abundances

One of the most important factors in the evolution of heavy-element abundances as a function of metallicity is the contribution from the s-process as well as the r-process. From the studies of solar composition at [Fe/H] = 0, more than 90% of strontium is considered to have been produced by the s-process. On the other hand, more than 90% of the europium in the solar system is believed to have been produced by the r-process. We have therefore investigated strontium and europium as typical elements that are produced mainly by the s-process or r-process, respectively.

The r-process component of the evolution of heavy-element abundances as a function of metallicity depends on how much each r-process site contributes. As seen in Section 3, the contribution of each r-process site depends on the element and also on the choice of the nucleosynthesis calculations and nuclear physics, where there are still large theoretical uncertainties. Therefore, we conducted a survey of the case where the relative contribution of each r-process site is artificially varied. For simplicity, we include only the collapsar and NSM contributions because similar results are obtained if we take CCSNe instead of collapsars in the metallicity region − 2 ≤ [Fe/H].

Figure 4 shows the evolution of the [Eu/Fe] ratio as a function of metallicity. The contribution fractions from NSMs are 100%, 90%, 50%, 10%, and 0% from left to right. The contribution from collapsars is the opposite. If less than 50% of europium is produced by NSMs, the evolution curve at [Fe/H] > −1 becomes flat as discussed in Côté et al. (2018). Although the contribution fraction of NSMs somewhat depends on the choice of the minimum coalescence times τ_g = 1, 10, and 100 Myr for binary NSMs, the flattening of the curve at [Fe/H] > −1 is independent of the choice of τ_g .

Zoom In Zoom Out Reset image size

As can be seen from Figures 1 and 2, in the presence of the NSM dynamical ejecta employing asymmetric FFD, the solar system's europium abundance is dominated by the NSMs. Since this is not consistent with the Galactic halo europium distribution, the fraction of the contribution from the dynamical ejecta component in NSMs is expected to be small.

On the other hand, it is very difficult to separate the contributions of collapsars and CCSNe in such a kind of schematic analysis. This is because the differences between these two components are seen in the ultra-metal-poor region, where the inhomogeneity is more pronounced. This is why the discussion in Section 3 is necessary.

Figure 5 shows the evolution of the [Sr/Fe] abundance ratio as a function of metallicity. The s-process contribution is taken from Kobayashi (2004). Again, from the analysis of the [Sr/Fe] distribution, a small NSM contribution is favored. However, this is not enough to distinguish which theoretical model illustrates the true contribution. Thus, a comparison between theoretical models and future observations is desired to provide better information as to which astrophysical sites a certain element comes from.

Zoom In Zoom Out Reset image size

To summarize, we have studied the cosmic evolution of the r-process abundance pattern in the context of GCE models that take into account multiple astrophysical sites simultaneously, i.e., NDW and MHDJ CCSNe, NSMs, and collapsars. The NSM r-process calculations were carried out with different input nuclear physics including symmetric and asymmetric FFDs. We then find that the r-process elements in the early galaxy are dominated by the yields from CCSNe and collapsars, while the NSM contribution is inevitably delayed due to the cosmologically long coalescence timescale for orbit decay by GW radiation. The relative NSM contribution rapidly increases with cosmic time. However, it does not reach even 1% of the total solar r-process composition until the metallicity is enriched to [Fe/H] ≥ −1.5. This conclusion does not change for a wide range of minimum coalescence times τ_g = 1–100 Myr in any GCE models including multiple sites and different input nuclear physics.

We also find that significant differences among our multiple-site GCE model calculations arise in the isotopic abundance pattern as a function of mass number A, while still satisfying the universality of elemental abundances for metal-poor halo stars. This is in contrast to previous studies that focused on only a single r-process site or a combination of at most two astrophysical sites in order to explain the universality of the elemental abundance pattern.

Several unique features of each astrophysical site are still expected in the GCE of the isotopic mass A abundance pattern. In particular, the collapsar contribution dominates from the very beginning of the early galaxy since its progenitor is a very massive star. The collapsar r-process shows an odd–even pattern over the entire mass range. Also, for both the collapsars and NSMs the abundance peaks shift toward the heavier mass region due to a high residual neutron-flux during the freeze-out of the r-process.

Although the elemental Z-abundance patterns are more or less similar to one another among the models, one can find exceptional differences in the actinides or light r-process elements, in particular Z < 42. Therefore, these are the important indicators of the dominant r-process site. Also, the peak structure is model dependent due to the different contributed fractions from the four astrophysical sites considered here.

It is therefore highly desirable to carry on spectroscopic observations with next-generation telescopes such as the Thirty Meter Telescope. These could provide the metallicity dependence of the isotopic abundance ratios of actinides, lanthanides, and lighter elements as well as the abundance peaks simultaneously. A separation of each r-process element into isotopes could provide constraints on the evolution of NSMs, MHDJs, and collapsars beyond that of models with a single r-process site.

This work is supported in part by Grants-in-Aid for Scientific Research of JSPS (20K03958, 17K05459). Work at the University of Notre Dame supported by DOE nuclear theory grant DE-FG02-95-ER40934. M.A.F. is supported by NASA grant 80NSSC20K0498.