{'fraction_non_alphanumeric': 0.04481260891050482, 'fraction_numerical': 0.0311061406018362, 'mean_word_length': 4.467861205915813, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 2, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Can neutron star mergers alone explain the r-process enrichment of the Milky Way?

Comparing Galactic chemical evolution models to the observed elemental abundances in the Milky Way, we show that neutron star mergers can be a leading r-process site only if at low metallicities such mergers have very short delay times and significant ejecta masses that are facilitated by the masses of the compact objects. Namely, black hole-neutron star mergers, depending on the black-hole spins, can play an important role in the early chemical enrichment of the Milky Way. We also show that none of the binary population synthesis models used in this paper, i.e., COMPAS, StarTrack, Brussels, ComBinE, and BPASS, can currently reproduce the elemental abundance observations. The predictions are problematic not only for neutron star mergers, but also for Type Ia supernovae, which may point to shortcomings in binary evolution models.


INTRODUCTION
Neutron star mergers (NSMs) 1 were proposed half a century ago (e.g., Lattimer & Schramm 1974) as a production site of the rapid neutron-capture process (r-process) elements such as europium, gold, and uranium. Recently, the gravitational-wave event GW170817 (Abbott et al. 2017a), also observed as a kilonova (AT2017gfo; Coulter et al. 2017;Smartt et al. 2017;Valenti et al. 2017;Tanvir et al. 2017;Abbott et al. 2017c) and a short gamma-ray burst (Abbott c.kobayashi@herts.ac.uk ilya.mandel@monash.edu et al. 2017b) with X-ray, optical and radio afterglows (Mooley et al. 2018;Lyman et al. 2018;Troja et al. 2017), confirmed NSMs as an r-process production site. The observational properties of AT2017gfo may be explained with neutron-rich dynamical ejecta, enhancing lanthanide production, and neutron-poorer outflows from the disc (Metzger & Fernández 2014;Cowperthwaite et al. 2017;Pian et al. 2017;Tanaka et al. 2017;Metzger 2019;Kawaguchi et al. 2018, but see Kasen et al. 2017), although the exact link between the kilonova observations and NSM ejecta components is not well understood. Strontium is directly detected in the day-1 spectra of the kilonova (Watson et al. 2019;Gillanders et al. 2022). Lanthanide production was supported by the day-2 spectra, but found to be rather small in later observational analyses (Waxman et al. 2018;Domoto et al. 2021). Con-straining the production of the third-peak elements such as gold and platinum is challenging.
Since then one more gravitational-wave detection of a neutron star (NS)-NS merger, GW190425, has been announced (Abbott et al. 2020a). In addition, two mergers of a NS with a black hole (BH) have been confidently detected, GW200105 and GW200115 (Abbott et al. 2021a); the provenance of GW190814, which could be a NS-BH or BH-BH merger, is unclear (Abbott et al. 2020b). The mass ejection from these BH-NS mergers is significant only for a small fraction of events (perhaps 10%, Drozda et al. 2020). There are simulations of BH-NS mergers (e.g., Janka et al. 1999;Rosswog 2005) but estimating the total mass ejection from BH-NS mergers requires 3D general relativistic (GR) hydrodynamical simulations with neutrino transport Siegel & Metzger 2018; see also Foucart et al. 2017 and the references therein). The mass of ejected material depends on the equation of state (EoS) of neutron stars, the spin of the BHs (Bauswein et al. 2014;Kyutoku et al. 2015) and the magnetic field (Paschalidis et al. 2015;Kiuchi et al. 2015), as well as the masses of the binary components (Krüger & Foucart 2020). The nucleosynthesis yields, however, have been provided only for a limited number of parameters in the numerical simulations (post-processing of hydrodynamical simulations; Bauswein et al. 2014;Wanajo et al. 2014;Just et al. 2015).
Galactic chemical evolution (GCE) models have been challenging NSMs as the main site of r-process nucleosynthesis (Argast et al. 2004;Cescutti et al. 2015;Wehmeyer et al. 2015;Haynes & Kobayashi 2019;Côté et al. 2019;Kobayashi et al. 2020a;van de Voort et al. 2020;Molero et al. 2021). The assumed delays between star formation and the occurrence of NSMs appear to be too long to explain the observed r-process enhancement (e.g., [Eu/Fe]) of extremely metal-poor stars in the Milky Way, while the rate of NSMs may be too low. The inclusion of inhomogeneous enrichment through state-of-the-art hydrodynamical simulations (Haynes & Kobayashi 2019;van de Voort et al. 2020van de Voort et al. , 2022 in some of these models leads to robust constraints on r-process enrichment, when the simulations predict the star formation and chemical enrichment history of our Milky Way Galaxy 2 ; inhomogeneous mixing does not help resolve discrepancies between model predictions and observations of stars. Therefore, galaxy modelling alone is unlikely to relieve the tension between observations and NSMs as the main r-process source within the current understanding of NSMs. Alternatively, r-process enrichment as-sociated with core-collapse supernovae has been developed: magnetorotational supernovae/hypernovae (e.g., Nishimura et al. 2015;Yong et al. 2021a;Reichert et al. 2023) or accretion disks/collapsars (Siegel et al. 2019, but see Just et al. 2022). Meanwhile, Grichener et al. (2022) suggested a new site in mergers of neutron stars with donor cores inside a common envelope.
The NSM rate has been estimated through binary population synthesis (BPS) models (e.g., Mennekens & Vanbeveren 2014Belczynski et al. 2016;Vigna-Gómez et al. 2018;Kruckow et al. 2018) 3 , whose predictions broadly match the merger rate inferred from Galactic binary pulsars, gravitational-wave events, and short gamma-ray bursts (see Mandel & Broekgaarden 2022 for a review). BPS models provide both a number of NSMs per unit star formation rate as a function of progenitor star metallicity and a delay-time distribution (DTD) between star formation and merger. The predicted DTD can be approximated with a simple t −1 power law, but there is significant variation between models and some metallicity dependence. The differences are due to a range of uncertain assumptions regarding all aspects of binary evolution: input parameter distributions, wind mass loss and stellar evolution, supernova remnant masses and natal kicks, the stability of mass transfer and the fraction of mass accreted by the companion, and the treatment of the common envelope phase (e.g., Broekgaarden et al. 2022).
The aim of this paper is not to perform a parameter study of r-process enrichment using GCE models and assuming NSM rates and DTDs as in previous works (e.g., Côté et al. 2019;Molero et al. 2021). Rather, we aim to derive an insight into binary evolution physics by considering what NSM yields and delay times would need to be in order to explain observed Milky Way r-process abundances with NSMs alone.
§2 briefly describes our GCE models. Model predictions and comparison to observational data are presented in §3. We discuss the implications for binary evolution models in §4. §5 gives our conclusions.

MODEL
We use the GCE code from Kobayashi et al. (2000) but including updated nucleosynthesis yields of AGB stars and core-collapse supernovae from Kobayashi et al. (2020a, hereafter K20) and those of Type Ia supernovae (SNe Ia) from Kobayashi et al. (2020b). The initial mass function (IMF) and the parameters for star formation histories are the same as in K20, which are determined to match other observational constraints in the solar neighborhood ( Fig. 1 of K20). In general these parameters are insensitive to elemental abundance tracks for constraining stellar physics (see Fig. A1 of Kobayashi et al. 2020b).
The core-collapse supernova yields are calculated based on stellar evolution and explosive nucleosynthesis by Umeda et al. (1999). For supernovae with explosion energies of ∼ 10 51 erg, a mass cut determines the remnant mass, while a mixing-fallback model is adopted for hypernovae (with energies > 10 51 erg), which mimics an aspherical explosion (Nomoto et al. 2013). These assumptions are chosen to match the resultant Fe yields with the observations of nearby supernovae (i.e., light curves and spectra; see §2.1.1 of K20), and the yield set can well reproduce the evolutionary trends of most elements.
The average evolution of chemical composition in the interstellar medium (ISM) is calculated by integrating the contributions of all relevant enrichment sources as a function of metallicity and timescale (see formulae in Kobayashi et al. 2000, or those of the single-star contribution in Matteucci 2021). The timescales are determined simply from the lifetimes of stars depending on their mass (and metallicity) for AGB stars and core-collapse supernovae, while the DTDs of merging binaries depend on many factors and model assumptions.
In this paper DTDs are defined as the rates of NSMs and SNe Ia per unit time per unit stellar mass formed in a population of stars with a coeval chemical composition and age (simple stellar population, SSP): R t = d 2 N/dt/dM * ,init , where M * ,init is the initial mass of the SSP.
There are some observational constraints on the DTDs in present-day galaxies (with observations biased toward high metallicity) or the rates integrated over time and volume containing various types of stars/galaxies (e.g., van den Heuvel & Lorimer 1996). However, for GCE, we need DTDs or rates as a function of metallicity of the SSP. BPS models (e.g., Tutukov & Yungelson 1993;Belczynski et al. 2002;Mennekens & Vanbeveren 2014;Kruckow et al. 2018;Vigna-Gómez et al. 2018) can provide the necessary information for GCE, i.e., the rate of events per unit stellar mass formed as a function of delay time and metallicity. However, predictions vary among the BPS codes because of uncertain assumptions regarding poorly understood binary physics, such as mass transfer including common envelope evolution or supernova natal kicks. De Donder & Vanbeveren (2004) were the first to explore such binary effects in a GCE code. Figure 1 shows the DTDs 4 of NSMs from various BPS codes, with the parameter sets recommended by the authors. In K20, for NSMs, the DTDs of the standard model from 4 Equivalently, the merger rates per unit stellar mass formed following a SSP. Mennekens & Vanbeveren (2014 for Z = 0.002 and Z = 0.02 were used, assuming a 100% binary fraction (which is also the case throughout this paper). Because of the larger mass of progenitors, BH-NS mergers have shorter minimal delay times than NS-NS mergers, and thus they are potentially good candidates for the early r-process enrichment. Moreover, at low metallicities, the rates of BH-NS mergers are predicted to be higher than those of NS-NS mergers, while at high metallicities predictions vary for different BPS simulations as follows. The BH-NS merger rate becomes lower than the NS-NS merger rate at high metallicities in Brussels, Kruckow et al. (2018, hereafter ComBinE), and Briel et al. (2022, hearafter BPASS), while the BH-NS rate is still slightly higher than the NS-NS merger rate in Mandel et al. (2021, hereafter COMPAS). For Belczynski et al. (2020, hereafter StarTrack), the BH-NS rate is so low and sparse that it is not plotted in this figure.
As in K20, the nucleosynthesis yields are taken from Wanajo et al. (2014, hereafter W14), and are calculated by post-processing a 3D-GR simulation of a NS-NS merger (1.3M + 1.3M ). Although the parameter dependence was discussed in K20 (see their §2.1.2), no yields were available for further investigation. In this paper, we also use nucleosynthesis yields from Just et al. (2015), which allow us to study the dependence of yields on NSM masses. While W14 yields included only dynamical ejecta of the NSM, Just et al. (2015) added the dominant BH torus ejecta driven by neutrinos and viscosity (α vis ) 5 . The spin of the BH formed after NSMs and the mass of the torus are also given as parameters. While the dynamical ejecta consistently have low electron fraction Y e and produce very heavy r-process elements, the secular (here called "wind") outflow component has relatively higher Y e and has more impact on element production with A < ∼ 130. For the EoS of NSs, we use the models 'SFHO' or 'DD2' if available in Just et al. (2015), and if not, we re-scale those with 'TM1' or 'TMA' by using the dynamical ejecta masses for 'SFHO' from Bauswein et al. (2013) and Bauswein et al. (2014).
The Eu yields are summarised in Table 1, compared with the value (labelled W14) used in K20. The ejecta mass is higher for systems with more extreme mass ratios, as shown in Fig. 7 of Bauswein et al. (2013). For equal-mass mergers, the mass of dynamical ejecta is similar between W14 and Just et al. (2015); it is 0.00483 M for the 1.35+1.35 M merger model in Bauswein et al. (2013) and 0.0143 M for 5 The ejecta from a remnant hypermassive NS, created in systems with lowenough binary mass and/or stiff-enough nuclear EoS, are not included. Depending on its lifetime and yet poorly understood magnetohydrodynamic processes of angular momentum transport and energy dissipation, the hypermassive NS could give rise to the ejection of substantial amounts of matter, presumably with relatively low neutron-richness compared to the other ejecta components (Martin et al. 2015;Fujibayashi et al. 2018).   (2015), depending on the viscosity parameter α vis ; the torus mass was also assumed to be 0.1, 0.03, or 0.3M respectively for the first three models, the massive NS-NS mergers, and the NS-BH mergers in Table 1 (see Just et al. 2015 for the other parameters of the remnant models). Although the results of Just et al. (2015) give ∼2-60 times larger Eu yields than W14, they also produce non-negligible amounts of Fe and O, and thus the [Eu/(O,Fe)] ratio is only larger by a factor of a few compared to W14. Hence, it is important to include not only the Eu yields but also the yields of other elements in GCE. When implementing these mass-dependent yields in GCE, we match BPS-modelled NSMs to the yields in Table 1 based on the total mass M tot ≡ M 1 + M 2 , the mass ratio q ≡ M 2 /M 1 , and the primary mass M 1 as described in the bottom section of the table. Whether or not BH-NS mergers eject significant mass and contribute to r-process nucleosynthesis depends on the BH mass, the mass and EoS of the NS, and on the spin of the BH and its tilt relative to the rotation axis of the binary (Bauswein et al. 2014;Foucart et al. 2018). If the BH mass is too high, NSs can directly plunge in without significant disruption; this was likely the case for both BH-NS mergers observed as gravitational-wave sources (Abbott et al. 2021a). We exclude BH-NS mergers with M 1 ≥ 6M from contributing to GCE 7 ; this limit holds for compact NSs, but a high BH spin and a stiff EoS of the NS could permit ejecta for even more massive BHs. The BH-NS models described in Table 1 assumed a BH dimensionless spin of 0.8; with a smaller spin, the mass ejection becomes much smaller. Following Bauswein et al. (2014), an additional parameter f spin is introduced to characterise the fraction of BH-NS mergers where the BH spin and tilt are favourable to sufficient matter ejection. We apply the nucleosynthesis yields (of all elements) of the BH-NS models in Table 1 to a fraction f spin of BH-NS mergers produced in BPS models in the relevant mass range.
When we use elemental abundances relative to Fe, the modelling of SNe Ia is also important. In K20, the progenitor model of Kobayashi & Nomoto (2009) is used, which is based on the theoretical calculation of optically thick winds from white dwarfs (WDs) by Hachisu et al. (2008). The DTDs depend on the metallicity of the binary systems, and are shown in Fig. 3 of Kobayashi & Nomoto (2009). These are similar to the "observed" DTDs from Maoz et al. (2014) at high metallicities (see also Fig. 12 of Kobayashi et al. 2020b). In this paper, we also use DTDs from BPS models for SNe Ia. Figure 1 shows the DTDs for single degenerate (SD) and double degenerate (DD) systems from various BPS codes, again with the parameter sets recommended by the authors. Even for a given BPS code, the parameter sets are not necessarily the same as for NSMs (see Di Stefano et al. in prep. for details). Note that Hachisu et al. (2008)'s WD wind effects are included only in Brussels. The metal-dependent nucleosynthesis yields of near-Chandrasekhar-mass (near-Ch-mass) and sub-Ch-mass models from Kobayashi et al. (2020b) are adopted for the SD and DD (or DD+He) systems, respectively. For ComBinE, DTDs of DD systems are given as a function of the total mass, and 1.37, 1.2, 1.1, 1.0, and 0.9M WD yields are used respectively for M tot > 2.12, 1.36, 1.16, 1.06, and 0.96M . For the other BPS, since the mass dependence is not available, we add 1.2, 1.1, 1.0, and 0.9M WD yields, respectively, with 10%, 40%, 40%, and 10% contributions toward the total sub-Ch-mass SNe Ia. The results are not very sensitive to this averaging.

RESULTS
Elemental abundances are being measured in the stellar atmospheres of nearby Milky Way stars for the past half century. Nowadays Galactic archaeology surveys with multiobject spectrographs are providing data for a million stars (e.g., Buder et al. 2021). However, in order to constrain stellar physics, elemental abundances with 3D non-localthermodynamic-equilibrium (NLTE) analysis of higher resolution data are necessary (Kobayashi et al. 2020a). Among r-process elements, Eu is the element accurately measured for the largest number of stars (see Fig. 32 of K20). Following Haynes & Kobayashi (2019), in order to remove the contribution from SNe Ia, we also show Eu abundances relative to α elements. Although Mg is the best observed α element, O is the element with the most robust nucleosynthesis yields (see §3.6 of K20) and is thus used in this paper. The black solid lines are the same as in Fig. 37 of K20, which can well reproduce the average evolution in all panels, thanks to magneto-rotational supernovae 8 , which produce most of Eu in this model. The other models consider NSMs as the only r-process site. Note that Brussels' DTDs for NSMs were already used in K20, but in this paper their DTDs of SD and DD SNe Ia are also used together with the metallicitydependent yields of Ch-mass and sub-Ch mass SNe Ia from Kobayashi et al. (2020b). Similarly, for the other BPS models, DTDs of both NSMs and all available channels of SNe Ia are used, without applying any arbitrary shift. Figure 2a shows that none of the BPS models can accurately reproduce the evolutionary change in [O/Fe] Fig. 2a.
Meanwhile, [Eu/O] evolution is shown in Figure 2c. The number of observational data points is reduced but the two best data sets remain, namely, the NLTE analysis from Zhao et al. (2016) Figure 2. This will help for some of the BPS models to match the solar ratios (at [Fe/H] = 0) but none of these BPS models can reproduce the observed evolutionary tracks; [Eu/(O,Fe)] should be increased more at early times than at present, if NSMs are the major r-process site.
In Figure 3, we include the BPS-predicted mass distribution of NSMs together with the mass-dependent yields (of all elements) from Just et al. (2015). We need to know the delay times (time from star formation to coalescence) and masses of the compact objects; we have such necessary information   for GCE only from COMPAS 9 and StarTrack 10 . The results depend on additional parameters: the Shakura-Sunyaev viscosity α vis of the BH torus plasma in Just et al. (2015) and fraction of BH-NS mergers that have spin and tilt suitable for sufficient mass ejection after mergers, f spin , in Bauswein 9 We use 2 models with Z = 0.0014 and 0.014 in  and exclude events with 'Optimistic CE'=True or 'Immediate RLOF>CE'=True. The COMPAS output already assumed a maximum neutron star mass of M NS,max = 2.0M . 10 We use the M30 model in Belczynski et al. (2020) but with Z = 0.0001, 0.0005, 0.001, 0.004, 0.008, 0.01, 0.02, and 0.03 and the events with 'hece'=0, which yields a NS-NS merger rate consistent with gravitationalwave observations. M NS,max = 2.5M is applied in the GCE code. et al. (2014). SNe Ia are also included with the same model as in K20. COMPAS has a high BH-NS merger rate and the nucleosynthesis predictions are thus highly sensitive to the f spin parameter. However, this model yields either too little Eu at low metallicities (with f spin = 0.1, blue lines) or too much Eu at high metallicities (with f spin = 0.5, cyan lines). α vis = 0.02 is assumed for NS-NS mergers in these two models.
StarTrack has a much lower BH-NS merger rate, and thus the result is insensitive to the spin+tilt parameter (f spin = 0.5 is assumed in these two models) but depends on the viscosity parameter. Higher α vis for NS-NS mergers (red lines) results in higher Eu abundances but these are still not large enough to match the observational data. The sharp increase in Eu for super-solar metallicities is also problematic. Note that the value of f spin estimated in BPS studies is much smaller (Drozda et al. 2020).
Having seen that existing BPS models struggle to reproduce the observations, we can return to our original question: If NSMs were the unique r-process site, what kind of DTDs could reproduce the observational data? To answer this question, we assume t −1 power-laws with fixed ranges for the DTD. The range and the height (i.e., rate) are arbitrarily determined in order to match the observations in Figure  4. Predictions for the following single power-law are shown with the red short-dashed lines: R NSM,t = 3.3 × 10 −7 t −1 M −1 for 0.01-20 Gyr R Ia,t = 3.33 × 10 −4 t −1 M −1 for 0.04-20 Gyr. (3.1) Here the 1.45+1.45 M yields are used, and the sub-Ch mass fraction to total SNe Ia is set to be 25%; these assumptions changes the normalization only. These single power-law functions cannot reproduce observations. The single powerlaw function gives a reasonable approximation for NSMs and SNe Ia from DD, but perhaps not for the other SNe Ia (Ruiter et al. 2011, for more discussion). The mismatches compared with the observations, apparent in Figure 4, are mainly caused by the minimum delay-times for both SNe Ia and NSMs; SNe Ia occur too early to keep the [α/Fe] plateau, and NSMs occur too late to create the [Eu/(O,Fe)] plateau. A flatter power than −1 gives a slightly shallower slope in the panels a and c, but it does not solve these problems. However, it is possible to reproduce the data if the DTDs strongly depend on metallicity (green long-dashed lines in Figure 4) as indicated by the gray solid lines in Figure 1:  As for Eq.3.1, the 1.45+1.45 M yields and the 25% sub-Ch mass fraction adopted. These functions given here are interpolated linearly in log Z at intermediate metallicities, while constant values are extrapolated at metallicities higher (lower) than Z = 0.02 (0.002).
For SNe Ia, optically thick winds from WDs can explain such a metallicity dependence of the rate (Kobayashi et al. 1998), although this effect is not fully included in existing BPS models. Metallicity dependent yields of Ch and sub-Ch mass SNe Ia (Kobayashi et al. 2020b) are also used throughout this paper.
For NSMs, nucleosynthesis yields at given NS/BH masses (with a given spin) are not expected to depend on the metallicity. The NSM DTD may strongly depend on the masses of compact objects, and one may think that the mass dependence could give a similar effect to the metallicity dependence. However, we find that the mass dependence alone (cyan dot-dashed lines) cannot reproduce the data in Figure  4. To show this, we assume an extreme case of DTD that strongly depends on the mass, using the 1.6+1.6 M and 1.4+5.08 M yields with f spin = 0.1, respectively, for NS-NS and BH-NS mergers, and 1.37M and 1.0M yields, respectively, for Ch and sub-Ch mass SNe Ia, with 25% sub-Ch mass fraction: Since BH-NS mergers yield more Eu than NS-NS mergers in Table 1, this model gives a higher [Eu/Fe] ratio at early times than the single power-law model. However, at later times, [Eu/O] keeps increasing and [Eu/Fe] shows a shallower decrease than observed. This is because star formation is ongoing, BH-NS mergers keep happening, and NS-NS mergers also start contributing at higher metallicities. Any changes of the range and normalization of the DTD do not help. It is not possible to reproduce the [O/Fe] ratios either; since our averaged Fe yield 11 from sub-Ch-mass SNe Ia is smaller than that of Ch-mass SNe Ia, the [O/Fe] ratio is higher than the single power-law model at [Fe/H]< 0, but is still lower than observed.
If BH spins strongly depend on metallicity, it becomes possible to reproduce the data. The blue dotted lines are for the mass-dependent DTDs in Eq.3.3, but assuming a higher f spin at low metallicities as follows: f spin = 0.1 for Z > 0.006 = 0.2 for Z ≤ 0.006.
(3.4) 11 The mass-dependent yields are averaged with the WD mass distribution in §2.
Here the fraction of sub-Ch mass SNe Ia of total SNe Ia is assumed to be only 5%; the coefficients are (9.5, 0.5), instead of (6, 2) in Eq.3.3. Note that the sub-Ch mass fraction changes the normalization only for the models with Eq.3.1, while it affects the shape of model curves with Eqs. 3.2 and 3.3. The sub-Ch mass fraction is very important for other elemental abundances, namely Mn, and found to be less than 25% from the Mn constraint in Kobayashi et al. (2020b).

DISCUSSION
The GCE models we use in this paper to explore the parameter space of NSMs can well predict the average evolution of elemental abundances but not the scatter. It is known that r-process elements show a large scatter particularly at low metallicities, which is caused by the rareness of the enrichment sources (e.g., Argast et al. 2004). There is also a selection bias of the observational data at [Fe/H] < ∼ − 3, where a few stars show very high r-process abundances, called r-II stars; the most metal-poor one (asterisk in Fig. 2) is the star in Yong et al. (2021a) at [Fe/H] = −3.5. Although these stars are well studied to obtain abundance patterns, the number of such stars is small (see Fig.14 of Yong et al. 2021b). Theoretically, the scatter can be predicted with stochastic GCE models (Cescutti et al. 2015;Wehmeyer et al. 2015), or more sophisticated, chemodynamical simulations of galaxies (Haynes & Kobayashi 2019;van de Voort et al. 2020van de Voort et al. , 2022 12 . Using the same yield set as in K20, Haynes & Kobayashi (2019) showed that the model with magneto-rotational supernovae reproduces the observed scatter of [Eu/(O,Fe)] well, while the model with NSMs shows a too large scatter with low [Eu/(O,Fe)] ratios. This is because the timescale of magneto-rotational supernovae is similar to that of Fe (and O)-producing core-collapse supernovae, while the NSMs are not associated with significant Fe (and O) production. Therefore when we compare the GCE models to the observed elemental abundances in Figures 2-4, it should be noted that the K20 model would scatter around the line, while the NSM-only models tend to scatter below the line.
Our analysis indicates that existing BPS models of NSMs struggle to reproduce the observed Eu abundances within our adopted yield models (Figs. 2 and 3). The metallicitydependent DTDs required to match observations are very different from those in existing BPS models (Figs. 1 and 4). Perhaps this should not come as a great surprise given that most BPS models do not accurately recover the Galactic double neutron star population properties. For example, the simulated mass distributions are often inconsistent with observed ones, which is likely a statement about imperfect recipes for translating the progenitor properties to remnant masses after supernovae (Fryer et al. 2012;Vigna-Gómez et al. 2018;. Meanwhile, Eu yields appear to depend very sensitively but uncertainly on the masses of merging objects, as shown in Table 1 and associated discussion, as well as on other imperfectly constrained quantities such as the electron fraction (Kullmann et al. 2022a), or on nuclear physics input such as fission fragment distribution and nuclear mass model (e.g., Eichler et al. 2015;Holmbeck et al. 2019;Kullmann et al. 2022b). Together, the likely flaws in population synthesis mass models and the uncertainty in NSM-mass-dependent yields make these comparisons challenging.
BPS models also struggle to reproduce the periodeccentricity distribution of Galactic double neutron stars, particularly the bifurcation in eccentricities for short orbital period systems (Andrews & Mandel 2019). This is, perhaps, more concerning than the mass distribution mismatch. Firstly, since the period and eccentricity of a double neutron star system after formation is directly linked to the delay time until merger through gravitational-wave emission (Peters 1964), this means that predicted DTDs are also likely to be inaccurate. Secondly and more fundamentally, it is an indication of significant uncertainties in the late stages of evolution preceding the formation of the second neutron star, particularly case BB mass transfer from a stripped expanding post He main sequence secondary onto the NS primary and the subsequent ultra-stripped supernova (Tauris et al. 2015), as well as an uncertainty in the natal kicks of NSs and lowmass BHs (e.g., Igoshev et al. 2021;Richards et al. 2022;Kapil et al. 2022).
In particular, the amount of mass stripped during case BB mass transfer, the dynamical stability of the mass transfer (whether or not the binary experiences common envelope evolution) and thus the ultimate binary separation, the remnant masses and the delay time until merger may well prove to depend on metallicity. The discovery of GW190425 (Abbott et al. 2020a), a NSM with a total mass significantly higher than any Galactic double neutron stars, hints at the possibility that there may exist a population of more massive NSMs with very tight orbits and short delay times and/or at lower metallicities, thus evading observation as Galactic radio pulsars but detectable as gravitational-wave sources (Safarzadeh et al. 2020;Galaudage et al. 2021;Vigna-Gómez et al. 2021; see also Kruckow 2020). For example, if lowermetallicity environments yield more massive stripped He stars after the first common-envelope phase, the subsequent case BB mass transfer would occur in a system with a more extreme mass ratio, leading to more binary hardening and a shorter delay time. Reduced mass loss in lower-metallicity environments could also produce tighter binaries and en-hanced tidal spin-up, leading to more rapidly spinning BHs (e.g., Bavera et al. 2021) in NS-BH mergers, supporting our metallicity-dependent f spin variation. However, many BPS models suggest that the BH in merging NS-BH binaries typically almost always forms when the binary is still wide (e.g., Broekgaarden et al. 2021;Chattopadhyay et al. 2022), reducing the opportunity for tidal spin-up of its progenitor.
This motivates the ad-hoc models proposed at the end of the previous section. The full space of allowed options is larger than we have considered: in addition to the freedom of choosing both mass-and metallicity-dependent DTDs, one could also vary the shape of the DTD (we did not explore this and restricted ourselves to t −1 power laws; see also Dvorkin et al. 2021). As mentioned above, more complex correlations, such as the possibility that low-metallicity environments yield both shorter DTDs and higher masses (and thus, perhaps, more Eu per NSM) are plausible. However, at least some of the models, such as the metallicity-dependent DTD model, are largely consistent with Eu abundance observations while retaining plausible parameters. For example, the total number of NSMs integrated over the age of the Universe is roughly 15-20 per million solar masses of star formation in both the single power law and metallicity-dependent DTD models. This is consistent with ∼ 30 NSMs per million years in the Galaxy, a value that matches inferred rates from Galactic radio pulsars (e.g., Pol et al. 2020), and with around 300 NSMs per Gpc 3 per year, which falls within the range of inferred rates from gravitational-wave observations (Abbott et al. 2021b).
The contribution of BH-NS binaries to Eu yields is even less certain. In order to contribute significant r-process nucleosynthesis, such mergers must disrupt the neutron star before it plunges through the innermost stable circular orbit (ISCO) of the BH, otherwise, insufficient material is left outside the BH. The radius of the ISCO scales linearly with the BH mass for a non-spinning (Schwarzschild) BH, but the tidal disruption radius scales only with the 1/3 power of the BH mass for sufficiently massive BHs. Therefore, if the BH is too massive, the neutron star will generally plunge in before it is disrupted. Disruption is easier if the BH is rapidly spinning and the inspiral is prograde, as the ISCO moves in (the ISCO radius can be up to a factor of 6 smaller for a maximally spinning Kerr BH than a Schwarzschild one). The condition on mass retention and nucleosynthesis is thus a function of the BH mass, spin and orbital tilt angle, and neutron star mass and EoS (Bauswein et al. 2014;Foucart et al. 2018). Given the uncertain mass and spin distribution of BHs in merging BH-NS systems (e.g., Broekgaarden et al. 2021), it is unclear whether a significant fraction of such mergers could contribute to r-process nucleosynthesis. Thus, while some short gamma-ray bursts and perhaps kilonovae have been conjectured to be associated with BH-NS merg-ers (Troja et al. 2008;Berger 2014;Li et al. 2017;Gompertz et al. 2020), other authors have argued that the contribution of BH-NS mergers to short gamma-ray bursts is likely low (e.g., Drozda et al. 2020). For example, it is probable that neither of the gravitational-wave observed BH-NS mergers would have experienced a kilonova (Abbott et al. 2021a).
Unfortunately, other observations are not yet sufficient to strongly constrain the merger masses and delay times of NSMs. Only a handful of mergers involving neutron stars have been so far observed with gravitational waves, making population inference challenging; furthermore, even inference on individual events may be sensitive to model assumptions (e.g., Mandel & Smith 2021). Kilonova observations not connected with gravitational-wave sources (e.g., Tanvir et al. 2013) could ultimately provide further constraints on both the NSM rates (see Mandel & Broekgaarden 2022 and references therein) and ejecta composition with more events, as expected in the Vera Rubin Observatory LSST data (Andreoni et al. 2019). Several attempts have been made to constrain DTDs using short gamma-ray bursts (e.g., Virgili et al. 2011;Wanderman & Piran 2015). However, such inference is sensitive to selection effects on both the prompt burst, X-ray afterglow and the later optical afterglow detection necessary for host galaxy identification. Short gamma-ray bursts can also be used to constrain the mass distribution of merging objects since a jet launch sets requirements on the masses, but these studies are not yet conclusive (Sarin et al. 2022;Salafia et al. 2022).
The presence of multiple uncertainties have led to explorations of the r-process abundances as constraints on the mass ratio distribution of NSMs (Bauswein et al. 2013;Holmbeck et al. 2021) or even on the NS EoS (Radice et al. 2018;Holmbeck et al. 2022). Placing stronger constraints will also require more detailed and faithful models of nucleosynthesis yields from mergers as a function of the masses of the merging objects. There are potentially sharp features in the ejecta masses and associated r-process abundances. The location of these features is not precisely known due to a combination of uncertainties in the merger models and the unknown EoS. Combined with the imprecise mass distribution predictions from BPS models, this creates an extra challenge for the kind of analysis undertaken here. Moreover, our notion of a "match" to observations that suffer from significant scatter is a somewhat imprecise visual averaging. This is partially due to the challenging observational selection effects when measuring abundances and partly to choices made in which stars to follow up. A homogeneous, non-biased observation of neutron capture elements for a wide range of metallicity (with NLTE analysis) would facilitate more robust comparisons between models and observations.

CONCLUSIONS
We investigate whether NSMs alone can explain the rprocess enrichment of the Milky Way. In particular, we explore whether NS-NS and especially NS-BH mergers in the early Universe can rapidly produce r-process elements instead of magneto-rotational supernovae. Provided that NSMs are the leading site of r-process nucleosynthesis, we constrain NSMs by comparing GCE models to the observed elemental abundances in the solar neighborhood, namely, As clearly indicated by the decreasing trend of [Eu/Fe] toward higher metallicities, DTDs of NSMs and SNe Ia should be very different. We show that none of the BPS used in this paper, i.e., COMPAS, StarTrack, Brussels, ComBinE, and BPASS, can reproduce the observations (Fig. 2); the DTDs of NSMs predicted from BPS have too long delay times, while those of SNe Ia have too short delay times.
We propose metallicity-dependent DTDs of SNe Ia and NSMs with t −1 power-laws that can explain the observations (Fig. 4). BH-NS mergers could play an essential role in the early chemical enrichment in the Milky Way if they dominate at low metallicities and are rare at high metallicities, and if a sufficient fraction of BH-NS mergers are able to produce significant mass ejections and contribute to r-process nucleosynthesis. However, the necessary BH-NS DTDs are also inconsistent with existing BPS models (Fig. 3).
CK thanks E. De Donder for discussions at a conference in Vulcano in 1999, which inspired some parts of this paper. CK acknowledges funding from the UK Science and Technology Facility Council (STFC) through grant ST/R000905/1 & ST/V000632/1. CK also acknowledges the Australian Research Council Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D), through project number CE170100013. The work was also funded by a Leverhulme Trust Research Project Grant on "Birth of Elements". IM acknowledges support from the Australian Research Council Centre of Excellence for Gravitational Wave Discovery (OzGrav), through project number CE17010004.