Total r-process Yields of Milky Way Neutron Star Mergers

While it is now known that double neutron star binary systems (DNSs) are copious producers of heavy elements, there remains much speculation about whether they are the sole or even principal site of rapid neutron-capture (r-process) nucleosynthesis, one of the primary ways in which heavy elements are produced. The occurrence rates, delay times, and galactic environments of DNSs hold sway over estimating their total contribution to the elemental abundances in the Solar system and the Galaxy. Furthermore, the expected elemental yield for DNSs may depend on the merger parameters themselves -- such as their stellar masses and radii -- which is not currently considered in many galactic chemical evolution models. Using the characteristics of the observed sample of DNSs in the Milky Way as a guide, we predict the expected nucleosynthetic yields that a population of DNSs would produce upon merger, and we compare that nucleosynthetic signature to the heavy-element abundance pattern of the Solar system elements. We find that with our current models, the present DNS population favors production of the lighter r-process elements, while underproducing the heaviest elements relative to the Solar system. This inconsistency could imply an additional site for the heaviest elements or a population of DNSs much different from that observed today.


INTRODUCTION
The trans-iron elements in the Solar system attribute the vast majority of their origins to neutron capture.Neutron capture is further divided into two main subclasses, slow ("s") and rapid ("r "), generally reflecting the neutron density of the environment in which each takes place (Burbidge et al. 1957;Cameron 1957).Requiring extreme neutron densities, the astrophysical site of the r -process has been a topic of debate for decades, with compact object mergers being proposed r -process sites since Lattimer & Schramm (1974).Since these early works, many studies have continued to claim that compact object mergers-particularly DNSs-can be responsible for the r -process (Lattimer & Schramm 1976;Symbalisty & Schramm 1982;Eichler et al. 1989;Freiburghaus et al. 1999).However, it was only recently that this theory has been observationally confirmed through the electromagnetic signature following the gravitational-wave detection of a binary neutronstar (NS) merger (NSM) (GW170817; Abbott et al. 2017a,b).Although the electromagnetic light curve and spectrum (AT2017gfo) associated with this event suggests that the r -process took place (Coulter et al. 2017;Cowperthwaite et al. 2017;Watson et al. 2019;Metzger 2019;Margutti & Chornock 2021), it is still unknown whether NSMs can account for the majority of the abundance of r -process elements observed in the Solar system and throughout the Galaxy (see Siegel 2022, for a review).Difficulties with the NSM model include explaining the r -process enrichment of ultra-faint dwarf galaxies with escape velocities far below the typical natal kick velocity applied to newborn NSs (Beniamini et al. 2016;Safarzadeh et al. 2019), the [Eu/Fe] abundances of metal-rich stars in the Milky Way disk (Côté et al. 2019), and inconsistencies between the observed lanthanide fraction of GW170817 and that of r -process enhanced stars in the Milky Way halo (Ji et al. 2019).Proposed additional r -process production sites include magnetorotational core-collapse supernova (Winteler et al. 2012;Nishimura et al. 2015), magnetized winds from proto-neutron stars (Thompson 2003), supernova fallback (Fryer et al. 2006), and accretion disk winds from collapsars (Surman et al. 2008;Siegel & Metzger 2018).
The study of NSMs in a Galactic context is a twopronged topic.First is the hydrodynamical investigation of the NSM itself and the detailed generation of r -process elements (Freiburghaus et al. 1999;Metzger et al. 2010;Roberts et al. 2011;Kasen et al. 2013;Rosswog et al. 2014;Grossman et al. 2014).Accurate simulations require a detailed treatment of hydrodynamics including general relativistic and radiative effects as well as a complex nuclear reaction network (for a review, see Metzger 2019).Second is the mixing processes in the interstellar medium and the broader chemical evolution of the Galaxy (Argast et al. 2004;Ishimaru et al. 2015;van de Voort et al. 2015;Côté et al. 2017;Macias & Ramirez-Ruiz 2018;Ojima et al. 2018;Macias & Ramirez-Ruiz 2019;Beniamini & Hotokezaka 2020;van de Voort et al. 2020;Wanajo et al. 2021).Studies disagree on the ability of turbulent mixing to efficiently bring r -process enriched material back to the Milky Way disk (e.g., Amend et al. 2022).However, these studies clearly show that the r -process site in the disk must be much rarer than normal core-collapse SN, consistent with NSM (Macias & Ramirez-Ruiz 2018, 2019), and that metal-rich stars such as the Sun must have been enriched by a large number of individual r -process events, so that stochastic effects can be ignored.
In this work, we assume the r -process abundance pattern of the Sun is representative of that of the Milky Way disk as a whole.Previously, Freiburghaus et al. (1999) and Metzger et al. (2010) have compared spectroscopic abundances of a variety of DNS merger simulations with the Solar abundance, finding reasonable agreement.However, the nucleosynthetic yield of NSMs depends on a combination of factors including the two components' masses, the NS Equation of State (EOS), and the various models for the contribution from disk winds and dynamical ejecta.While the latter three may be uncertain, the known population of DNSs in the Milky Way inform us about the masses of merging DNSs.Here, we combine a variety of merger simulations with a realistic population of DNSs informed by the observed Milky Way systems, accounting for those systems that have already merged and are not seen to-day, to determine a realistic nucleosynthetic yield of rprocess elements at the present day.This study differs from previous studies of Galactic chemical evolution in that we explore NSM yields that depend on the NS masses and radii.We also systematically explore the effect of nuclear uncertainties on our NSM yields.We begin our study by calculating the r -process abundances that could be expected from the Milky Way population of DNSs, starting with an overview of the observed sample of DNSs and its nuances in Section 2. We then describe our procedure for calculating nucleosynthetic yields from a DNS merger with a particular combination of NS masses in Section 3. In Section 4 we combine the DNS sample with our model of elemental yields to derive a nucleosynthetic fingerprint from a Galaxy-scale population of NSMs.Finally, we provide some conclusions in Section 5.

DNS SYSTEMS
In the decades since the detection of the first DNS, the Hulse-Taylor binary (Hulse & Taylor 1975), the number of observed systems has steadily increased.Today, there are ≃20 known or suspected DNSs in the Milky Way (for a review, see Tauris et al. 2017).Of these, roughly half are expected to merge due to gravitational wavedriven orbital decay within a Hubble time.Throughout this work, we ignore DNSs with orbits too wide to merge within a Hubble time since there is evidence that these systems formed through a different evolutionary subchannel (Andrews & Mandel 2019).After further restricting for only systems with well-measured masses, we arrive at a sample of seven DNS1 which we list in Table 1 along with their masses.
Not all of the seven DNSs in our sample are equally weighted, which can be explained through two broadly defined effects (Narayan et al. 1991;Phinney 1991;Kalogera et al. 2001Kalogera et al. , 2004)).First, the shorter-period systems merge quickly, such that for every system observed, many more have already formed and merged.A complete analysis accounts for this effect by weighting each system by the inverse of the DNS lifetime, estimated to be the sum of the pulsar characteristic age and the time to merger.(This weighting scheme also makes apparent why systems that take more than a Hubble time to merge have a negligible contribution to the characteristics of the merging DNS population.)Second, observational biases need to be taken into account.For example, intrinsically fainter systems with a smaller Table 1.Observed DNS systems and their ejecta properties for our fiducial model.The median and 16th and 84th percentiles between all random samples is shown for the wind and dynamical ejecta masses (mi) and their respective mass fractions of lanthanides and actinides (X lan i ).3.9 +1.9 −1.9

Name
1.9 detectable volume contribute more significantly, as each detected system implies a larger underlying population.
Properly incorporating observability requires a detailed treatment of pulsar lifetimes, beaming fractions, pulsar luminosities, and the footprints and depths of pulsar surveys within which they were discovered.A quantitative analysis combining all of these effects for the seven DNSs in our sample was undertaken by Pol et al. (2019Pol et al. ( , 2020)).Although additional systems have been detected since these works (e.g., Sengar et al. 2022), none will merge within a Hubble time.
In Table 1 we include the relative weighting of each DNS which is based on the derived merger rate reported by Pol et al. (2019).
Four of the seven systems (J0737−3039, J1906+0746, B1913+16, and J1757−1854) have a significantly higher weight, due principally to the larger differences in gravitational wave-driven merger times.
We additionally consider B2127+11 (also known as M15C), which is found within the core-collapsed globular cluster M15 and therefore likely to be formed dynamically.Although this system will merge in a Hubble time, dynamically formed systems are not expected to contribute significantly to the overall merger rate of DNSs in the universe; decades of focused radio observations of globular clusters implies limited numbers of such systems, in agreement with theoretical models (Ye et al. 2019(Ye et al. , 2020)).We treat this system separately in our subsequent analysis, removing it from any extrapo-lation to the r -process abundance pattern of the Milky Way.
Finally, we include the two DNSs detected by the LIGO-Virgo network of gravitational wave observatories, GW170817 (Abbott et al. 2017a) and GW190425 (Abbott et al. 2020).Analysis of the gravitational waveforms leads to mass measurement accuracies of ∼0.1 M ⊙ .In our subsequent analysis, the two stars' masses (the mass ratio in particular) dictate the resultant abundance yield.While GW170817 has component masses consistent with the Milky Way field population, GW190425 is unique.Its masses are sufficiently discrepant from any other known DNS leading to multiple theories about its formation (Romero-Shaw et al. 2020;Safarzadeh et al. 2020;Vigna-Gómez et al. 2021).We consider these gravitational wave detections as well, but leave them out of any derivation of the predicted Milky Way-averaged r -process abundance pattern.

NSM YIELDS
We consider two main sources of NSM ejecta: dynamical and disk-wind.The total elemental yield from a merger will be a product of the two ejecta masses and their respective compositions.We employ a forward method in which both the ejecta masses and compositions are dependent on the NS masses and the EOS.We combine the contributions from both the dynamical ejecta and the disk wind into a total yield for each of the 10 observed DNS systems listed in Table 1.First we describe our fiducial model in Section 3.1, then we de-scribe how we approximate uncertainties on these models in Section 3.2.We discuss the effects of alternatives to our fiducial prescriptions in Sections 4.2 and 4.3.For a schematic overview of our model, see Figure 3 in Holmbeck et al. (2022).

Fiducial Model
For our fiducial model, we adopt the DD2 EOS (Hempel & Schaffner-Bielich 2010).To describe the dynamical ejecta mass, we use Equation 6of Krüger & Foucart (2020), which is a fit to a series of numerical relativity simulations.We opt for the prescription provided by Dietrich et al. (2020) in Equation S4 for the mass of the accretion disk as a function of NS masses and the EOS, which is based on a fit to a series of 73 separate numerical relativity simulations performed by different groups.For details on why we choose these prescriptions as our default choices, see Holmbeck et al. (2022).We will explore the uncertainty associated with these choices in Section 4.2.
Next, the masses of each of these components are multiplied by their respective calculated elemental composition.We use the nuclear reaction network code Portable Routines for Integrated nucleoSynthesis Modeling (PRISM; Sprouse 2020) with the Finite Range Droplet Model (FRDM2012 Möller et al. 2012) nuclear mass model to calculate the nucleosynthetic yield of the NSM ejecta.
For the dynamical ejecta we use a single trajectory from NSMs from the 1.4-1.4M ⊙ NSM simulations of S. Rosswog (Piran et al. 2013;Rosswog et al. 2013) as in Korobkin et al. (2012) with varying initial electron fraction, Y e (see the calculations in Holmbeck et al. 2019Holmbeck et al. , 2021)).Rather than use a single Y e to describe the composition of the dynamical ejecta, we allow the mass to be distributed in a narrow Gaussian around a central Y e value.This initial central Y e is chosen from the average Y e estimate from Equation 10 in Nedora et al. (2021), which is written as a function of NS mass ratio and tidal deformability.
For the disk wind component, the elemental yield is affected by the lifetime of the NSM remnant, both in the composition and the total amount of mass that is ejected.We take the remnant lifetime to depend on the merger parameters as in Lucca & Sagunski (2020) and represent the fraction of the disk mass that will be ejected by fitting data in Metzger & Fernández (2014) as follows: where τ is the NSM remnant lifetime.Then, we use the lifetime-based nucleosynthesis yields and fitting formulae as in Holmbeck et al. (2021) to obtain a lifetimedependent composition.In this way, a longer-lived merger remnant will eject more mass and be generally less neutron-rich than its short-lived counterpart.
Once both the component ejecta masses (M ej ) and yields (Y ) are calculated, they are simply multiplied together to obtain a total abundance yield.

Model Uncertainties and Error Propagation
Significant uncertainty surrounds each step in our forward model.We have chosen reasonable prescriptions for the dynamical and wind ejecta masses and compositions, but they are certainly not the only plausible descriptions for what are intrinsically complex, multiscale physical processes.We therefore estimate the uncertainty in our model taking a Monte Carlo approach; when calculating the expected yields for each NSM, we take 50,000 random samples overall to account for the effects of multiple separate sources of uncertainty.
First, for the Galactic and globular cluster DNS pairs, we take 50 random draws for the two NS masses from a normal distribution with centroids and standard deviations given by the values in Table 1.Because the NS masses in the Galactic DNS systems are measured with high precision, the systematic uncertainty on the elemental yields from these measurements is a negligible effect.However, the NS masses for the for the gravitational wave-detected systems are a significant source of uncertainty.For the two systems GW170817 and GW190425, instead of sampling from a normal distribution, we take 50 samples from their posterior mass distributions: the low-spin prior mass pairs from Abbott et al. (2019) 2 and Abbott et al. (2020)3 , respectively.
For each DNS, we put the 50 NS mass pairs through the separate equations that determine the dynamical and wind ejecta masses, the remnant lifetimes, and the dynamical ejecta centroid Y e .Next, 1000 samples are drawn randomly from normal distributions that describe the uncertainty associated with each of these four quantities.Specifically, the uncertainty on the ejecta mass can be >50% and have significant effect on the calculated yields.Uncertainty on the abundances depends on the uncertainty in τ and Y e for the wind and dynamical ejecta, respectively.For the disk lifetime (and therefore the disk ejecta mass), the median uncertainty is 100%.The centroid of the dynamical ejecta Y e is highly sensitive to the NS masses, significantly affecting .Predicted r -process yields of merger ejecta for DNS systems in the Galaxy compared to the r -process abundance pattern of the solar system, scaled by a constant (gray).Colored lines show individual random samples, while the black lines show the average between all draws.Note that B2127+11 is a system found in a globular cluster (see text for details).the dynamical ejecta yields.Therefore, for each DNS in Table 1, we propagate and quantify our systematic uncertainty by calculating 50,000 sets of yields.

RESULTS
For each DNS system modeled with our fiducial model, Table 1 shows the median wind and dynamical ejecta masses, their respective lanthanide mass fractions, and the 16 th and 84 th percentile for each of these quantities, represented as upper and lower uncertainties.Overall, these DNSs eject roughly equal amounts of matter owing to their similar masses.Notably, the wind ejecta mass dominates over the dynamical ejecta mass by about a factor of 50.In addition, note that the lanthanide mass fraction of the wind ejecta is 1.9×10 −3 for almost all mass pairs.This phenomenon can be attributed to the effect of the relatively stiff EOS on the remnant lifetime.The remnant lifetime prescription adopted in our fiducial model (Lucca & Sagunski 2020) is highly sensitive to the maximum non-rotating NS mass, M TOV .The somewhat high M TOV of the DD2 EOS (≈2.4 M ⊙ ) leads to lifetimes consistently above 530 ms, causing the NSM to eject its entire disk and result in lower neutron richness for the r -process in the wind ejecta (see, e.g., Lippuner et al. 2017 andHolmbeck et al. 2021).Note, our fiducial model also therefore predicts that GW170817 was longlived (τ > 530 ms), in agreement with more rigorous studies that find lifetimes for the GW170817 remnant O(1 s) (e.g., Gill et al. 2019).GW190425 is altogether different.The significantly larger NS masses and higher mass asymmetry cause the NSM to eject very little material through the wind (three orders of magnitude lower than that of GW170817), though we note the wind has a large lanthanide mass fraction in our model.
In addition to the overall yields of the individual ejecta components, our method allows us to calculate projected detailed elemental abundance patterns.Figure 1 shows the expected yields for each merging Galactic DNS system, as well as the DNS found in the globular cluster M15, B2127+11.Each colored line shows one of the ran-dom samples drawn within the model uncertainties.For plotting simplicity, we only show 200 randomly selected samples from the full set of 50,000 samples.The black lines display the element-by-element averages calculated from the entire set of random draws.The abundances are scaled by the total ejecta mass (see Section 3.1), but note that they are not yet scaled by the weight.Interestingly, the yields are directly comparable, and it can be directly seen that the systems ejecta roughly equal amounts of r -process material.However, this total yield is difficult to directly compare to observational abundances due to the unknown mass in which this yield may be diluted.Therefore, we also plot the Solar system r -process pattern (gray line) for reference, which we vertically shift by the same, arbitrary value for each system.
The low total masses and high symmetry values of the Galactic DNSs typically lead to strong production of the first r -process peak (e.g., 38 Sr) relative to the rest of the heavy elements.J1913+1102, the only DNS system in this sample with significantly unequal neutron star masses, produces the highest amount of actinides and the highest total ejecta mass overall.The predicted abundances for the globular cluster system, B2127+11, behave similarly to the majority of DNSs in the Milky Way, though it has one of the lowest yields of actinide abundances.
Predicted abundance patterns for the DNSs detected by gravitational waves from their merger are shown in Figure 2. Recall that for these systems, the neutron star masses are drawn as co-dependent pairs from the (low-spin) posterior mass distributions, as described in Section 3.This spread on the masses is an additional source of uncertainty on the abundance patterns and predicted ejecta masses, leading to a wider spread in the abundance patterns that are produced.Interestingly, of all the systems, GW170817 shows the highest actinide abundance.In addition, its total ejecta mass is among the largest of the DNSs in this sample, which- Figure 3.Total ejecta mass and lanthanide mass fractions for each DNS in Table 1.
like J1913+1102-can be attributed to its slightly higher mass and asymmetry than DNSs in the Galaxy.In contrast to all other systems we consider here, we find that the abundance pattern for GW190425 is suppressed in the first r -process-peak elements.Additionally, the high NS masses and mass ratio of GW190425 is predicted to eject very little mass, though with very high lanthanide mass fraction, possibly due in part to its low first r -process-peak mass.Compounded with its relatively large distance and poor sky localization, the low luminosity produced from this small amount of ejecta mass could in part explain why a kilonova counterpart could not be found (Antier et al. 2020;Coughlin et al. 2019;Hosseinzadeh et al. 2019;Lundquist et al. 2019).
The combined total ejecta masses and lanthanide mass fractions for the DNSs are displayed in Figure 3.For most of the DNSs, the total ejecta masses are consistently large (≳0.1 M ⊙ ), while the lanthanide mass fractions are ≲10 −2 .These DNSs are consistent with GW170817, in which our model predicts a total lanthanide mass fraction of ≃ 10 −2.25 , in agreement with several model studies of the AT2017gfo kilonova (Arcavi et al. 2017;Chornock et al. 2017;Kilpatrick et al. 2017;McCully et al. 2017;Nicholl et al. 2017;Tanaka et al. 2017;Tanvir et al. 2017;Troja et al. 2017).On the other hand, Figure 3 highlights how our fiducial model predicts vastly different ejecta masses (≲ 10 −2 M ⊙ ) and lanthanide fractions (≳0.1) for GW190425.

Fiducial model: total yields
With the individual yields for each Milky Way DNS calculated, we can also test the combined yields from these systems.For this combination, we take the masses from the seven Galactic DNS systems and weigh their contribution by how representative each DNS is of similar-mass NSMs in population synthesis models.The weights are listed in Table 1, and details of their calculation can be found in Andrews (2020).We multiply each system's total yield by their respective weight in Table 1 and divide by the sum of the total weights.As described in Andrews ( 2020), the weight is essentially the expected merger rate (in Myr −1 ) for each type of system.
Figure 4 shows this combined, weighted abundance patterns of the seven Galactic DNSs.Although lowmass, very symmetric DNS systems are the most common in the Milky Way sample, higher-mass systems like B1913+16 are expected to merge more frequently and thus contribute more to the total abundance.As in Figures 1 and 2 we compare the resulting abundance pattern to that of the Solar system in gray.We again emphasize that the Solar system r -process pattern is representative of the relative pattern of each of element; the absolute vertical scaling is arbitrary, as there is no accurate way to convert the Solar system abundance pattern to an individual merger ejecta pattern.Combined, the total merger yields from Milky Way DNS systems look similar to the Solar system r -process pattern, though with a higher enhancement to the lighter of the r -process elements.

Macroscopic model uncertainties
We have chosen state-of-the-art nuclear reaction and decay rates, calculated self-consistently within the FRDM2012 nuclear model as well as some of the most up-to-date fits to hydrodynamic models of the disk and dynamical ejecta.In addition, the composition of the disk ejecta was found by calculating the nucleosynthesis from tens of thousands of individual tracer particles from hydrodynamic simulations that vary the lifetime of the merger remnant.Nevertheless, there are several additional areas of systematic uncertainty in this analysis.For example, the dynamics of the accretion disk wind remain a topic of intense study with few to no constraints.Furthermore, other analytic descriptions of the dynamical and post-merger ejecta masses exist (e.g., Dietrich & Ujevic 2017;Radice et al. 2018;Coughlin et al. 2019;Dietrich et al. 2020;Krüger & Foucart 2020;Nedora et al. 2021).We test 25 unique combinations of dynamical and disk ejecta masses from the literature and found that there was very little sensitivity on the total combined abundance pattern.In all model ejecta variations we tested, the long lifetime of the merger remnant caused the wind component to consistently dominate over the dynamical ejecta.The largest variations we observed were that the dynamical ejecta mass described by Radice et al. (2018), which leads to a slightly less-abundant third r -process peak, and the description by Dietrich & Ujevic (2017), which produces a higher third peak.However, we find the uncertainty from the disk and dynamical ejecta descriptions to be minor overall, especially compared with the microscopic model assumptions which we discuss in the subsequent section.

Microscopic model uncertainties
While the nucleosynthetic yield is relatively insensitive to the exact choice of disk and dynamical ejecta models for our fiducial model, these variations become important when changing the EOS and nuclear mass model.To demonstrate these effects, we choose a softer EOS, SFHo (Steiner et al. 2013), and show the combined Milky Way abundances from the seven Galactic DNSs in the left panel of Figure 5 and report the results in Table 2.The lower M TOV and smaller NS radii given with SFHo leads to much smaller wind ejecta masses (ranging from a factor of five to 500) and dynamical ejecta masses consistently about two times larger then the results in Table 1.The lanthanide mass fractions also vary for the components.Generally with the SFHo model, the wind ejecta can reach higher neutron richness (allowed by lower remnant lifetimes), with a maximum lanthanide mass fraction of 10 −1.32 (B1913+16).On the other hand, the average Y e for the dynamical ejecta is predicted to be somewhat higher with SFHo, causing the dynamical ejecta lanthanide mass fraction to be overall comparable to, but slightly lower than, our fiducial EOS choice (DD2).The effect of the softer EOS is that first and second r -process-peak production is somewhat decreased, leading to a more consistent overall fit with the Solar r -process abundances.We summarize the combined ejecta masses and lanthanide mass fractions with the SFHo model in Figure 6.Note that there is more variation in the ejecta masses and lanthanide production when using SFHo, unlike in Figure 3 in which the calculations tended to group around the same values for nine of the ten DNS systems.
Next, we run all of our nucleosynthesis calculations with the Duflo-Zuker nuclear mass model (DZ; Duflo & Zuker 1995) instead of our baseline FRDM2012 model.For this test, we revert back to the DD2 EOS.Note that changing the nuclear mass model does not affect the total NSM ejecta mass.This change required recalculating thousands of thermodynamic trajectories, which we are able to accomplish using the PRISM code.The right panel of Figure 5, shows the independent effect of changing the nuclear mass model on the combined abundance pattern.The net effect of employing a different nuclear mass model for the composition of the ejecta is that the lanthanides are produced in higher amounts and the actinides are slightly reduced.In addition, the r -process peaks have lower abundances relative to the baseline case in which FRDM2012 nuclear masses are used.

DISCUSSION AND CONCLUSIONS
In this work, we have presented our model predictions for the r -process abundances from NSMs inside and outside of the Galaxy.We find that the total contribution by Galactic NSM systems produce overall lower actinide abundances compared to gravitational-wave detected NSMs (i.e., compare Figures 4 and 2).In addition, our predictions for r -process nucleosynthesis from Milky Way DNSs fail to match the Solar actinide abundances, which may imply that DNSs are not sufficiently frequent or high-yield to account for what we see in the Galaxy on average.We have demonstrated that this conclusion is robust to reasonable changes to choices for the macroscopic model, i.e., changes to the disk wind and dynamical eject models.On the other hand, the microscopic model uncertainties produce significant changes to the resulting elemental yields.When using a softer EOS, SFHo, we find an elemental abundance pattern generated by NSMs more closely matching the Solar pattern, but a lanthanide fraction for GW170817 too high compared with observations.Similarly, employing a different nuclear mass model (DZ) for the nucleosynthesis calculations produced overall higher yields of the first and second r -process peaks.
While we have not exhaustively examined the entire range of possible model variations, none of our models can simultaneously fit the entire Solar r -process abundance pattern.With either a different EOS describing the ejecta or a different nuclear mass model affecting the composition, the actinides are consistently underproduced relative to the Solar system lanthanides.This inconsistency suggests that a contribution from events that produce higher actinide abundances may be needed: either more asymmetric NSMs, a contribution from neutron-star black-hole mergers, or other sites that were active in the Galaxy's early history.
In addition to the multiple model uncertainties we adopt, our approach intrinsically makes several assumptions in extrapolating to the entire r -process production history of the Milky-Way.First, when we combine our individual NMS yields into one "Galactic DNS" abundance pattern, we inherently assume that the ejecta has mixed homogeneously into the interstellar medium (ISM) and in proportion to both the total ejecta mass and the weight ascribed to each DNS.Fully capturing the subtleties of NSM rates and yields and how those quantities interact with star-formation demands a full galactic chemical evolution analysis, which is beyond the scope of this work.Our combined NSM yields in Figure 4 can instead be thought of as a representation of the average Milky Way r -process contribution by NSMs rather than a robust prediction.
Second, the Galactic DNSs may not be a representative sample of a population of merging neutron-star binaries.The pulsars in our sample of Galactic DNSs in Table 1 all have relatively young spin-down ages (Tauris et al. 2017), suggesting they were made from metalrich stars.Massive metal-poor binaries, like those seeding the Milky Way with r -process elements early in its formation, may produce different mass neutron stars and therefore a different elemental abundance pattern upon merger.Additionally, the detection of GW190425 confirms that neutron star binaries more massive than those currently observed in the Milky Way sample can exist.Although its formation is yet to be uncovered, two possibilities have been presented in the literature to explain the mass discrepancy between this system and other Galactic DNSs; either a formation channel produces DNSs at very short orbital periods that they quickly merge, making them unlikely to be observed (e.g., Romero-Shaw et al. 2020), or some process (such as mass fallback, burying their magnetic fields) exists that keeps the NSs in some systems from becoming pulsars, or at least forms pulsars with a short lifetime (Safarzadeh et al. 2020;Vigna-Gómez et al. 2021).Either possibility suggests that the observed sample listed in Table 1 forms an incomplete, or "hidden," representation of the population of merging DNSs.Nevertheless even if a hidden population exists, Figure 3 shows that massive NSMs such as GW190425 produces at least an order of magnitude less ejecta mass compared with the other systems.Assuming that the rate of NSMs from any putative hidden population does not dwarf the merger rate of those systems observed in the Milky Way, we can be optimistic that any subpopulation of GW190425-like systems missing from the current Milky Way census would provide at most a perturbation to the integrated Milky Way abundance pattern.Finally, any comparison between the model data and the Solar values inherently assumes that NSMs are the main source of elements from the first r -process-peak to the actinides.We have neglected additional sources from, e.g., exotic SNe that may produce a robust rprocess pattern.These sites are thought to be roughly ten times more frequent than NSMs, but may eject ten times less the material that an NSM would (e.g., Tsu- Figure1.Predicted r -process yields of merger ejecta for DNS systems in the Galaxy compared to the r -process abundance pattern of the solar system, scaled by a constant (gray).Colored lines show individual random samples, while the black lines show the average between all draws.Note that B2127+11 is a system found in a globular cluster (see text for details).

Figure 2 .
Figure 2. Predicted abundance patterns of merger ejecta for NSMs observed by LIGO/Virgo.

Figure 4 .
Figure 4. Predicted combined abundance patterns of merger ejecta for DNSs in the Galaxy.

Figure 5 .
Figure 5.Effect on the combined abundance pattern of merger ejecta for DNSs in the Galaxy from a different EOS (left) and a different nuclear mass model (right).

Figure 6 .
Figure6.Total ejecta mass and lanthanide mass fractions using the SFHo EOS.Note that GW190425 is off the scale of this plot near Mej ∼ 10 −4 and X lan ∼ 10 −1 .

Table 2 .
Observed DNS systems and their ejecta properties as in Table1, but for the SFHo EOS.