A new delay time distribution for merging neutron stars tested against Galactic and cosmic data

TL;DR

This paper introduces a new delay time distribution for merging neutron stars, tested against gamma-ray burst data and Galactic chemical evolution, providing insights into their role in element production and cosmic events.

Contribution

The paper derives a novel delay time distribution for neutron star mergers and tests it against multiple observational constraints, improving understanding of their timing and contribution to element synthesis.

Findings

SGRB distribution fits with 300-500 Myr delay times

Eu production timescale must be <30 Myr if MNS are sole Eu sources

Metallicity dependence and core-collapse SNe contribute to Eu production

Abstract

The merging of two neutron stars (MNS) is thought to be the source of short gamma-ray bursts (SGRB) and gravitational wave transients, as well as the main production site of r-process elements like Eu. We have derived a new delay time distribution (DTD) for MNS systems from theoretical considerations and we have tested it against (i) the SGRB redshift distribution and (ii) the Galactic evolution of Eu and Fe, in particular the [Eu/Fe] vs [Fe/H] relation. For comparison, we also tested other DTDs, as proposed in the literature. To address the first item, we have convolved the DTD with the cosmic star formation rate, while for the second we have employed a detailed chemical evolution model of the Milky Way. We have also varied the DTD of Type Ia SNe (the main Fe producers), the contribution to Eu production from core-collapse SNe, as well as explored the effect of a dependence on the metallicity of the occurrence probability of MNS. Our main results can be summarized as follows: (i) the SGRB redshift distribution can be fitted using DTDs for MNS that produce average timescales of 300-500 Myr; (ii) if the MNS are the sole producers of the Galactic Eu and the occurrence probability of MNS is constant the Eu production timescale must be on the order of <30 Myr; (iii) allowing for the Eu production in core-collapse SNe, or adopting a metallicity-dependent occurrence probability, allow us to reproduce both observational constraints, but many uncertainties are still present in both assumptions.

Tables

Table 1. Table 1: The table reports the value of the occurrence probability α M ​ N ​ S subscript 𝛼 𝑀 𝑁 𝑆 \alpha_{M\!N\!S} for each tested DTD that must be inserted in order to obtain the present time cosmic rate of MNS at z = 0 𝑧 0 z=0 equal to the gravitational wave event rate found by the LIGO/Virgo Collaboration (Abbott et al., 2017a ) . In the third and fourth columns are reported the fractions of MNS systems merged after 20 Myr and 100 Myr from a single initial starburst.

MNS	${\alpha }_{MNS}$	merged	merged
DTD	($\times {10}^{-2}$)	in 20 Myr	in 100 Myr
Our DTD $\beta =0.9$	$0.66$	0.008	0.056
Our DTD $\beta =0.0$	$0.81$	0.019	0.129
Our DTD $\beta =-0.9$	$1.02$	0.039	0.241
Our DTD $\beta =-1.5$	$1.18$	0.056	0.336
Prop to $t^{-1}$	$1.09$	0.096	0.319
Prop to $t^{-1.5}$	$1.60$	0.301	0.703
Prop to $t^{-2}$	$1.73$	0.500	0.901
Constant 10 Myr	$1.75$	1.000	1.000

Table 2. Table 2: The table reports all the models shown in the graphs. As explained in section 5, the occurrence probabilities for SNeIa ( α I ​ a subscript 𝛼 𝐼 𝑎 \alpha_{Ia} ) has been tuned to obtain the current Galactic rate of SNeIa and the current cosmic rate MNS. Instead, the occurrence probability α M ​ N ​ S subscript 𝛼 𝑀 𝑁 𝑆 \alpha_{M\!N\!S} for MNS has been taken from table 1. The yield of Eu from MNS ( 7 t ​ h superscript 7 𝑡 ℎ 7^{th} column) has been tuned to reproduce the abundance in our Sun (the results of our simulations are reported in the 9 t ​ h superscript 9 𝑡 ℎ 9^{th} and the 10 t ​ h superscript 10 𝑡 ℎ 10^{th} columns). Finally, in the 8 t ​ h superscript 8 𝑡 ℎ 8^{th} column is reported the predicted current rate of MNS in the Milky Way.

Name	Eu from	SNeIa	MNS	${\alpha }_{Ia}$	${\alpha }_{MNS}$	MNS Eu yield	Current MNS	$X_{Fe}$	$X_{Eu}$
	CC-SNe	DTD	DTD	($\times {10}^{-2}$)	($\times {10}^{-2}$)	($\times {10}^{-6}{M}_{\odot }$)	rate (${\text{Myr}}^{-1}$)	($\times {10}^{-3}$)	($\times {10}^{-10}$)
1A	no	Prop to $t^{-1}$	Our DTD $\beta =0.9$	4.77	0.66	4.0	76	1.38	1.68
2A	no	Prop to $t^{-1}$	Our DTD $\beta =0.0$	4.77	0.81	4.0	84	1.38	2.29
3A	no	Prop to $t^{-1}$	Our DTD $\beta =-0.9$	4.77	1.02	4.0	94	1.38	3.15
4A	no	Prop to $t^{-1}$	Our DTD $\beta =-1.5$	4.77	1.18	4.0	102	1.38	3.77
5A	no	Prop to $t^{-1}$	Prop to $t^{-1}$	4.77	1.09	4.0	93	1.38	3.27
6A	no	Prop to $t^{-1}$	Prop to $t^{-1.5}$	4.77	1.60	4.0	104	1.38	4.65
7A	no	Prop to $t^{-1}$	Prop to $t^{-2}$	4.77	1.73	4.0	93	1.38	4.28
8A	no	Prop to $t^{-1}$	Constant 10 Myr	4.77	1.75	4.0	94	1.38	4.29
1B	no	G05 ${\beta }_a=-0.9$	Our DTD $\beta =0.9$	4.67	0.66	4.0	76	1.52	1.68
2B	no	G05 ${\beta }_a=-0.9$	Our DTD $\beta =0.0$	4.67	0.81	4.0	84	1.52	2.29
3B	no	G05 ${\beta }_a=-0.9$	Our DTD $\beta =-0.9$	4.67	1.02	4.0	94	1.52	3.15
4B	no	G05 ${\beta }_a=-0.9$	Our DTD $\beta =-1.5$	4.67	1.18	4.0	102	1.52	3.77
5B	no	G05 ${\beta }_a=-0.9$	Prop to $t^{-1}$	4.67	1.09	4.0	93	1.52	3.27
6B	no	G05 ${\beta }_a=-0.9$	Prop to $t^{-1.5}$	4.67	1.60	4.0	104	1.52	4.65
7B	no	G05 ${\beta }_a=-0.9$	Prop to $t^{-2}$	4.67	1.73	4.0	93	1.52	4.28
8B	no	G05 ${\beta }_a=-0.9$	Constant 10 Myr	4.67	1.75	4.0	94	1.52	4.29
1AS	yes	Prop to $t^{-1}$	Our DTD $\beta =0.9$	4.77	0.66	1.5	76	1.38	2.57
2AS	yes	Prop to $t^{-1}$	Our DTD $\beta =0.0$	4.77	0.81	1.5	84	1.38	2.80
3AS	yes	Prop to $t^{-1}$	Our DTD $\beta =-0.9$	4.77	1.02	1.5	94	1.38	3.13
4AS	yes	Prop to $t^{-1}$	Our DTD $\beta =-1.5$	4.77	1.18	1.5	102	1.38	3.22
5AS	yes	Prop to $t^{-1}$	Prop to $t^{-1}$	4.77	1.09	1.5	93	1.38	3.11
6AS	yes	Prop to $t^{-1}$	Prop to $t^{-1.5}$	4.77	1.60	1.5	104	1.38	3.68
7AS	yes	Prop to $t^{-1}$	Prop to $t^{-2}$	4.77	1.73	1.5	93	1.38	3.55
8AS	yes	Prop to $t^{-1}$	Constant 10 Myr	4.77	1.75	1.5	94	1.38	3.56
1BS	yes	G05 ${\beta }_a=-0.9$	Our DTD $\beta =0.9$	4.67	0.66	1.5	76	1.52	2.57
2BS	yes	G05 ${\beta }_a=-0.9$	Our DTD $\beta =0.0$	4.67	0.81	1.5	84	1.52	2.80
3BS	yes	G05 ${\beta }_a=-0.9$	Our DTD $\beta =-0.9$	4.67	1.02	1.5	94	1.52	3.13
4BS	yes	G05 ${\beta }_a=-0.9$	Our DTD $\beta =-1.5$	4.67	1.18	1.5	102	1.52	3.22
5BS	yes	G05 ${\beta }_a=-0.9$	Prop to $t^{-1}$	4.67	1.09	1.5	93	1.52	3.11
6BS	yes	G05 ${\beta }_a=-0.9$	Prop to $t^{-1.5}$	4.67	1.60	1.5	104	1.52	3.68
7BS	yes	G05 ${\beta }_a=-0.9$	Prop to $t^{-2}$	4.67	1.73	1.5	93	1.52	3.55
8BS	yes	G05 ${\beta }_a=-0.9$	Constant 10 Myr	4.67	1.75	1.5	94	1.52	3.56

Equations

$$R_{MNS}(t)=k_{\alpha }{\int }_{{\tau }_i}^{\text{min}(t,{\tau }_x)}{\alpha }_{MNS}(\tau )\psi (t-\tau )f_{MNS}(\tau )d\tau$$

$${\int }_{{\tau }_i}^{{\tau }_x}{f}_{MNS}(\tau )d\tau =1.$$

$$k_{\alpha }={\int }_{m_L}^{m_U}\varphi (m)dm$$

$${\tau }_{gw}=\frac{0.15A^4}{m_1{m}_2(m_1+m_2)}Gyr$$

$$\mu (M_{DN},m_2)=M_{DN}^2\cdot m_2-m_2^2\cdot M_{DN},$$

$$\mu (M_{DN})=0.25M_{DN}^3$$

$${\tau }_{gw}=\frac{0.6A^4}{M_{DN}^3}Gyr.$$

$$df({\tau }_{gw})=df(A,M_{DN})=g(A)h(M_{DN})dAd{M}_{DN},$$

$$\{\begin{array}{l}g(A)\propto A^{\beta }\\ h(M_{DN})\propto const.\end{array}$$

$$f(\tau_{gw})d\tau_{gw}=d\tau_{gw}\int_{A_{l}}^{A_{s}}g(A)\,h(M_{D\!N})\,\Bigl{|}\frac{\partial M_{D\!N}}{\partial\tau_{gw}}\Bigr{|}\,dA$$

$$A_{x}=\Bigl{(}\frac{M_{D\!N,x}^{3}\tau_{gw}}{0.6}\Bigr{)}^{1/4}$$

$$f({\tau }_{gw})\propto \frac{1}{{\tau }_{gw}^{4/3}}\frac{1}{\beta +7/3}[A^{\beta +7/3}{]}_{A_l}^{A_s}.$$

$$f({\tau }_{gw})\propto {\tau }_{gw}^{(1/4)\beta -3/4}(M_{DN,s}^{(3/4)(\beta +7/3)}-M_{DN,l}^{(3/4)(\beta +7/3)}).$$

$$f(\tau )\propto \{\begin{array}{l}0\text{if}\tau <10Myr\\ p_1\text{if}10<\tau <40Myr\\ p_2{\tau }^{0.25\beta -0.75}(M_{DN,s}^{0.75(\beta +2.33)}-M_{DN,l}^{0.75(\beta +2.33)})\\ \text{if}40Myr<\tau <13.7Gyr\end{array}$$

$$\Psi (z)=\frac{1+p_{1}z}{1+(z/z_p{)}^{p_2}}$$

$$\Psi (z)\propto \{\begin{array}{l}(1+z{)}^{-3.08}z<1.60\\ (1+z{)}^{-4.98}z\ge 1.60\end{array}$$

$$\Psi (z)=\frac{0.015(1+z{)}^{2.7}}{1+((1+z)/2.9{)}^{5.6}}{M}_{\odot }{\text{Mpc}}^{-3}{\text{yr}}^{-1}$$

$$\varphi (m)=0.171m^{-2.35},$$

$$\begin{array}{rl}\frac{d{G}_i}{dt}(r,t)& =-\psi (r,t)Z_i(r,t)+\\ & +{\int }_{m(t)}^{M_U}\psi (t-\tau (m))Q_{mi}(t-\tau (m))\varphi (m)dm+\\ & +Z_{i,0}A(r,t)\end{array}$$

$$\log (Fe/H{)}_{\odot }=-2.792\log (Eu/Fe{)}_{\odot }=-6.496$$

$$X_{Eu}=3.50\times 10^{-10}{X}_{Fe}=1.34\times 10^{-3}.$$

$${\alpha }_{MNS}(Z)=\{\begin{array}{ll}3{\tilde{\alpha }}_{MNS}& Z_{Fe}\le -1.0\\ {\tilde{\alpha }}_{MNS}& Z_{Fe}>-1.0\end{array}$$

$${\alpha }_{MNS}(Z)=\{\begin{array}{ll}{\tilde{\alpha }}_{MNS}(1+5\ln (-Z))& Z_{Fe}\le -1.0\\ {\tilde{\alpha }}_{MNS}& Z_{Fe}>-1.0\end{array}$$

Full text

††pubyear: 2019††pagerange: A new delay time distribution for merging neutron stars tested against Galactic and cosmic data–A new delay time distribution for merging neutron stars tested against Galactic and cosmic data

1 Introduction

The merging of two neutron stars via the emission of gravitational waves has been invoked to explain a variety of phenomena like the production of heavy r-process elements (Korobkin et al., 2012; Hotokezaka et al., 2013) and the short gamma-ray bursts (SGRBs) (Eichler et al., 1989; Narayan et al., 1992). The observation of binary systems composed by compact and massive star remnants (Hulse & Taylor, 1975; Tauris et al., 2017) with measurable changes in their orbital periods, and the recent detection of the gravitational wave (GW) transient GW170817 (Abbott et al., 2017a) strongly support the occurrence of this process.

Elements heavier than Fe cannot be produced by exoenergetic fusion reactions in stars and instead they must be the result of neutron capture on Fe-peak nuclei. If the ingestion of neutrons is "rapid" compared to the $\beta$ decay timescales, heavy nuclei are produced, while "slow" captures result in relatively lighter products (Burbidge et al., 1957; Seeger et al., 1965). Correspondingly, these elements are called r- and s-process elements. Moreover, r-process elements can be subdivided in a light (with A $<$ 90) and a heavy (A $>$ 90) subclass. A widespread used tracker for heavy r-process elements is Europium (Eu) (e.g. Matteucci et al., 2014). The production of heavy r-process elements requires very neutron rich environments, such as core-collapse SNe (CC-SNe) and MNS. The dominant production site, though, remains uncertain (Côté et al., 2018a).

Early studies, based on the observation of r-process elements in very metal poor ([Fe/H] $<$ -3.0 dex) stars, led to the conclusion that r-process production should occur in massive, short-lived stars (m $>$ 10 $M_{\odot}$) (Truran, 1981).

Hydrodynamical simulations (e.g. Woosley et al., 1994; Wheeler et al., 1998) supported such a result, showing that the neutrino-driven wind from newly born neutron stars (NS) in stars with masses higher than 20 $M_{\odot}$ is a promising (although imperfect) site for r-process nucleosynthesis. However, subsequent studies (Arcones et al., 2007; Arcones & Janka, 2011; Martínez-Pinedo et al., 2012) have cast doubts on the capability of this mechanism to produce a robust r-process abundance pattern. In particular, it was found that neutrino-driven winds were not enough neutron rich to produce elements with A $>$ 120. Moreover, the final abundance pattern is very sensitive to details of the physical conditions at explosion. A rare sub-class of metal poor CC-SNe, called magneto-rotational (MR) SNe, has been theorized to be able to overcome these problems (Cameron, 2003; Winteler et al., 2012), though their ability to reproduce the solar abundances of these elements is not clear (Woosley & Heger, 2006), but see also (Cescutti & Chiappini, 2014).

In parallel, other authors have explored the role of compact object mergers (CBM) (like NS-NS and NS-BH systems) in the production of r-process elements. Nucleosynthesis calculations (Korobkin et al., 2012; Eichler et al., 2015) predict the ejection of up to $10^{-2}\,M_{\odot}$ of r-process matter with a robust abundance pattern in a single coalescence event, due to the large number of neutrons per target Fe nuclei. Therefore, this source appears to be more reliable with respect to the CC-SNe.

The viability of double neutron star mergers as r-process production site has been tested in galactic chemical evolution simulations, which trace the abundances of different elements in the gas of the Milky Way. In this regard, a very important parameter to evaluate is the timescale of production of a given element, which can be estimated as the average coalescence timescale of a single burst stellar population. For MNS systems, the delay time (or time between the birth of the binary system and its final merging) is determined by both the stellar nuclear lifetime and the gravitational delay time, and it can vary in a wide range. The typical timescale for enrichment from MNS depends on the distribution of the delay times (DTD).

A popular choice in literature consist in selecting a DTD in the form of a power law, i.e. $\propto t^{\gamma}$ with $0.5\leq\gamma\leq 2.0$ (Côté et al., 2017b; Hotokezaka et al., 2018) where $t$ is the total (i.e. nuclear+gravitational) coalescence delay time; in other studies, the total coalescence delay time is obtained by adding a constant gravitational delay on top of stellar evolutionary lifetimes (Matteucci et al., 2014; Vangioni et al., 2015). Some authors, like Matteucci et al. (2014), have tested both MNS only, and a combination of MNS and CC-SNe as Eu producers. In all studies it is found that, in order to reproduce the decreasing trend in the [Eu/Fe] vs [Fe/H] relation in the thin disk stars, the Europium pollution timescale should not exceed 100 Myr or be even shorter (10-30 Myr).These timescales can be attained either using very steep power laws ($\propto t^{-2}$) or assuming a constant gravitational delay of 1 Myr on top of the nuclear evolution timescale of the progenitors. These DTDs are in disagreement with the results of population synthesis models (Belczynski et al., 2017; Giacobbo & Mapelli, 2018a; Côté et al., 2018a) that predict a $\gamma$ in the range $[-1.0;-1.5]$.

Gamma-ray bursts (GRBs) are luminous transients of high energy photons (in the range of 10-100 keV) incoming isotropically to the Earth at a rate of 1 per day (Klebesadel et al., 1973). They can be divided in short ($<2$ s, SGRBs) and long ($>2$ s, LGRBs) (Norris et al., 1984; Meszaros, 2003); SGRBs are characterized by high peak energy, while LGRBs exhibit a softer emission. It is generally accepted that short and long GRBs are produced by different mechanisms; in particular the latter have been explained as due to the acceleration of matter along the poles of a newly born black hole during massive CC-SNe explosions (MacFadyen & Woosley, 1999). The SGRBs are instead thought to be produced from the magnetic recombination or neutrino-antineutrino annihilation during the merging of compact binary objects (Narayan et al., 1992).

Despite the small statistics, various groups (Ghirlanda et al., 2016; Zhang & Wang, 2018) have reconstructed the redshift distribution of SGRBs, while others (D’Avanzo, 2015; Fong et al., 2017) have tried to reproduce it by convolving different DTDs with the cosmic star formation history. It has been found that the SGRB redshift distribution is better reproduced when a DTD $\propto t^{-1}$ or $\propto t^{-1.5}$ is used. Moreover, between 1/4 and 1/3 of the total SGRBs comes from early-type galaxies (Berger, 2014; Fong et al., 2017), which host old stellar populations. This evidence contrasts with the notion of short timescales to describe MNS, in favor of longer ones.

A great deal of information comes from the detection of GW170817. Not only it has been the first, direct evidence of a MNS (Abbott et al., 2017a), but it was also associated with a SGRB (Abbott et al., 2017b) and an optical transient known as "kilonova" (Smartt et al., 2017). The study of its multiband lightcurves and of the spectra has provided constraints on the composition of the ejecta, though a consensus has not yet been reached. Smartt et al. (2017) concluded that these ejecta should have been rich in light r-process and very deficient in heavy r-process, while others (Tanvir et al., 2017; Troja et al., 2017; Evans et al., 2017; Pian et al., 2017) support the hypothesis of a robust production of extremely opaque r-process elements. Finally, based on this event, Abbott et al. (2017a) estimated the cosmic local rate of MNS events to be $1540_{-1220}^{+3200}$ $\text{yr}^{-1}\text{Gpc}^{-3}$.

In this paper we want to test under which conditions it is possible to accomodate the various observational constraints under the hypothesis that MNS are responsible for all SGRBs, and that they are important contributors to the Europium chemical enrichment of the Milky Way. To do this we: (i) derive a new formulation for the DTD from theoretical considerations, adapting the Greggio (2005) procedure, which was developed for Type Ia Supernovae, to the progenitors of MNS; (ii) use the cosmic star formation history (CSFR) and the redshift distribution of SGRBs to evaluate the fraction of massive stars which give rise to a MNS event; (iii) test the results versus the Milky Way abundance pattern of Eu and Fe using a full galactic chemical evolution model of the Milky Way, which successfully reproduces a majority of observational constraints (e.g. Matteucci et al., 2009). We also explore the effect of adopting the DTD formulations which are proposed in the literature.

2 Our new DTD

The delay time distribution function $f_{M\!N\!S}(\tau)$ for MNS represents the distribution of the coalescence times from an instantaneous burst of star formation of unitary mass and it is fundamental to compute the MNS rate. The MNS rate can be written as:

[TABLE]

where $\psi$ is the SFR, $k_{\alpha}$ is the number of stars with mass in a suitable range per unit mass of star forming gas in a stellar generation, and $\alpha_{M\!N\!S}$ is the fraction of them which gives rise to MNS events. In principle, both $k_{\alpha}$ and $\alpha_{M\!N\!S}$ can vary as a function of time, but for the sake of simplicity we will assume them constant. In Sect. 5.4 we will explore the effect of varying $\alpha_{M\!N\!S}$. The time $\tau$ is the delay time defined in the range $(\tau_{i},\tau_{x})$, where $\tau_{i}$ is the minimum delay time of occurrence for merging neutron star (here fixed to 10 Myr) and $\tau_{x}$ is the maximum delay time which can be larger than a Hubble time, depending on the progenitor model. Since $f_{M\!N\!S}(\tau)$ is a distribution function, it must be normalized to 1 in the allowed range for $\tau$:

[TABLE]

The parameter $k_{\alpha}$ depends on the IMF, and it is:

[TABLE]

where $m_{L}$ and $m_{U}$ are, respectively, the progenitor minimum and maximum mass to produce a NS. In this paper we adopt $m_{L}=9M_{\odot}$ and $m_{U}$ = 50 $M_{\odot}$.

We now derive the form of $f_{M\!N\!S}(\tau)$. For MNS the time delay is the sum of the evolutionary lifetime of the secondary component of the binary system plus the gravitational time delay. From Landau et al. (1962):

[TABLE]

where $m_{1}$ and $m_{2}$ are the masses of the primary (more massive) and of the secondary star of the binary system, respectively, while $A$ is the initial separation of the neutron star binary system. We can now show that this equation can be rewritten as a function of the total mass, instead of the masses of the two components, with a small error. The mass dependent term, namely the denominator of (4), can be written as:

[TABLE]

with $M_{D\!N}=m_{1}+m_{2}$ i.e. the total mass of the system. This function can be studied in its two variables. For a fixed $M_{D\!N}$ this formula produces a parabola with the maximum at $m_{2}=0.5M_{D\!N}$.

The lightest known neutron star whose mass has been precisely measured is the companion of the binary pulsar PSR J0453+1559, with an estimated mass of $1.17\pm 0.004$ $M_{\odot}$ (Martinez et al., 2015). On the other hand, the heaviest measured neutron star is the PSR J0348+0432, with a mass of $2.01\pm 0.04$ $M_{\odot}$ (Antoniadis et al., 2013). In a yet not fully confirmed study, Linares et al. (2018) estimate the mass of PSR J2215+5135 to be $2.27\pm 0.15$ $M_{\odot}$, and this fact would rule out some of the proposed equations of state for NS interior (Özel & Freire, 2016).

We have decided to adopt a slightly different range with a lower limit to the mass of the neutron star of 1 $M_{\odot}$ and an upper limit of 2 $M_{\odot}$. Therefore the allowed range for $M_{D\!N}$ goes from 2.0 to 4.0 $M_{\odot}$. The requirements on the mass of the primary NS further constrains the range of possible values for $m_{2}$ so that not all combinations of $M_{D\!N}$ and $m_{2}$ are acceptable. In figure 1 we show the function $\mu(M_{D\!N},m_{2})$ plotted vs $m_{2}$ for five different values of $M_{D\!N}$. The shaded are shows the portion of the plane in which $m_{1}>m_{2}$. The sections of the parabolae within the shaded area are very flat: for each $M_{D\!N}$ the function $\mu$ varies very little (no more than 10%) inside the allowed region of the parameter space ($m_{1}$,$m_{2}$).

We can conclude that it is possible to rewrite $\mu$ as a function of $M_{D\!N}$ with a negligible error: every binary system with the same total mass and the same initial separation has approximately the same delay. Substituting the value of the maximum for each $M_{D\!N}$, $\mu(M_{D\!N},m_{2})$ becomes:

[TABLE]

and so it follows that:

[TABLE]

Having simplified the Landau formula we derive distribution of the delay times due to the emission of gravitational waves as follows.

The contribution to the number of systems with delay $\tau_{gw}$ from systems with separation $A$ and total mass $M_{D\!N}$ can be written as:

[TABLE]

where $g(A)$ is a function describing the initial separations, while $h(M_{D\!N})$ is the distribution of the total masses of the systems that will merge. In Eq. (8) these two functions are assumed to be independent. We further assume that $g(A)$ can be described as a power law and, for simplicity, that the distribution of binary masses is flat. Therefore the two functions can be written as:

[TABLE]

As can be seen, $\beta$ parametrizes the shape of the distribution in initial separations. Integrating Eq. (8) over the relevant range of separations we derive the number of MNS systems which merge with a gravitational time delay $\tau_{gw}$:

[TABLE]

with $A_{l}$ and $A_{s}$ are the minimum and maximum separations that contribute to $n(\tau_{g\!w})$ for systems with total mass $M_{D\!N}$. They can be calculated as:

[TABLE]

where the subscript x stands for s and l. The result of such calculation is:

[TABLE]

Substituting Eq. (11) in Eq. (12) we derive the distribution of the gravitational wave delay times as:

[TABLE]

The number of systems with a total delay between $\tau$ and $\tau+d\tau$ will be close to the number of systems with a gravitational delay between $\tau_{gw}$ and $\tau_{gw}+d\tau_{gw}$, because of the short and similar nuclear lifetime ($\tau_{n}$) of stars with high masses, that ranges between $\sim 30$ Myr for 9 $M_{\odot}$ stars and $\sim 4$ Myr for 50 $M_{\odot}$ stars. However, the existence of a distribution of $\tau_{n}$ implies that the DTD deviates from $f(\tau_{gw})$ at the short delay times, because the parameter space is limited from the requirement $\tau=\tau_{n}+\tau_{gw}$. Greggio (2005) developed a DTD for binary white dwarfs taking this effect into account. Following her results we construct the following distribution of total (nuclear + gravitational) delay time:

[TABLE]

where $p_{1}$ and $p_{2}$ will be chosen so as to obtain a continuous and normalized function. The first portion of the distribution ends with the formation of the first double neutron star system. Ten million years is the nuclear lifetime of a typical massive star. The second portion refers to systems which merge soon after the formation of the double neutron star binary system. Similar to the distributions in Greggio (2005), this portion of the distribution function is described as a flat plateau, up to the maximum nuclear time delay of the double neutron star binary or, in other words, the lifetime of the minimum mass progenitor of a neutron star ($\sim 9$ $M_{\odot}$). The third portion is just the distribution of the gravitational delay times and attains to those systems for which the time delay is dominated by the gravitational radiation mechanism described by the Landau equation. The results can be seen in figure 2 where a normalization factor of 1 has been assumed.

2.1 Other tested DTDs

As a comparison for DTDs, we have also tested four other common options found in literature. Three of them are single-slope power laws, with $\gamma=-1.0,\,-1.5,\,-2.0$, as proposed by Côté et al. (2018b, a). We again impose a minimum value to the total delay time of 10 Myr, which represents the shortest lifetime of neutron stars progenitors. All of them are normalized between this effective minimum time and one Hubble time (13.7 Gyr). The various choices of DTDs discussed so far are shown in figure 2.

We also tested the DTD implemented in Matteucci et al. (2014), in which the total delay time is equal to 1 Myr on top of the stellar lifetime, implicitely assuming that all binary neutron stars are born very close, and merge soon after they are formed. In Matteucci et al. (2014) it has been shown that longer gravitational delays do not reproduce the Eu abundances in very metal poor stars.

2.2 The progenitors of SNeIa

Since one of our constraints concerns the trend of [Eu/Fe] vs [Fe/H] and since an important part of the iron enrichment comes from SNeIa, we need to specify a choice for their DTD. We test two different options:

1. [i.]

a distribution proportional to $t^{-1}$, as suggested by Totani et al. (2008) and Maoz et al. (2012) on the basis of empirical data; 3. 2.

the DTD in the wide double degenerate scenario depicted in Greggio (2005) with an initial separation distribution coefficient $\beta_{a}=-0.9$ and a maximum nuclear lifetime $\tau_{n,x}=0.4$ Gyr. The reasons for these choices are presented in Matteucci et al. (2009).

Both DTDs are normalized in the range 0.04-13.7 Gyr, and we show them in figure 3. These two distributions are not very different. In fact, at times larger than $\sim 3$ Gyr, the fraction of exploded SNeIa is similar. On the other hand, as will be show in section 5, their effect on the [Eu/Fe] vs [Fe/H] relation at low (-2.5 < [Fe/H] < -1.0 dex) metallicities is pronounced. Chemical evolution models based on DTDs very different from this kind of shape have difficulties to fit the abundances of low metallicity stars.

The various DTDs discussed so far imply different timescales for the Europium and Iron pollution on the ISM. As mentioned in the Introduction, the pollution timescale can be represented as the average coalescence timescale for a single burst stellar population, or the time within which half of the events occur. Figure 4 shows the cumulative distribution of the DTDs for the MNS and the SNeIa’s adopted in our computations. The Europium pollution timescales for the various models are very different, ranging from $\sim 20$ Myr to 2 Gyr as the DTD for the MNS changes from the steepest pure power law to the flattest of our proposed DTD. Instead, in our models, the timescale for Fe production from SNeIa ranges only from 0.4 to 0.7 Gyr.

3 Comparison with the SGRB redshift distribution

As briefly described in the Introduction, it has been possible to derive the SGRB redshift distribution despite the small statistics. However, there are discrepancies between the results of different groups.

3.1 The SGRB redshift distribution

To constrain our models we consider two redshift distributions of SGRB proposed in the literature: the Ghirlanda et al. (2016) (herefater G16) and Zhang & Wang (2018) (hereafter ZW18).

G16 considered seven observables: (i) peak flux, (ii) fluence, (iii) observer frame duration, (iv) observer frame peak energy distribution, (v) redshift, (vi) isotropic energy and (vii) isotropic luminosity, usually not available at the same time for each event, and a set of parametric relations that bind these together. The parametric form chosen for the redshift is:

[TABLE]

Then, they generated synthetic populations of SGRBs, choosing randomly the values of the parameters, and compared their features to the observational dataset. A total of 10 free parameters (3 for redshift distribution, 3 for peak energy distribution and 4 for energy-luminosity correlations) have been tested. This is referred in the paper as ’model a’. G16 considered also a slightly different scheme for relations between observables with independent (from the peak energy and between themselves) distributions for luminosity and duration, raising the number of free parameters to 11. This is referred as ’model c’.

The values of the redshift distribution parameters which best represent the observational characteristics of the SGRB population turn out to be: $p_{1}$ = 2.8 (model a) or 3.1 (model c), $p_{2}$ = 3.5 (model a) or 3.6 (model c) and $z_{p}$ = 2.3 (model a) or 2.5 (model c).

With a different approach, ZW18 first considered the 16 SGRBs with known redshift. From them, they derived a peak energy - luminosity relation and have used this relation to find the redshift of the other 284 SGRBs in the Fermi-GBM sample. Finally, ZW18 have used these inferred redshifts to fit a bimodal relation in the form of:

[TABLE]

This is an ever-decreasing function of redshift, sharply in contrast with the result of G16.

3.2 Simulations and results

In this section we compare the SGRB redshift distribution shown above to theoretical expectations of the redshift distribution of MNS based on the DTDs described in Sect. 2. To do this we convolve the DTDs with the cosmic star formation rate (CSFR) (see Eq. (1)), for which we adopt the relation obtained by Madau & Dickinson (2014):

[TABLE]

base on a Salpeter (1955) IMF, that is:

[TABLE]

when normalized between 0.1 and 100 $M_{\odot}$.

For a Salpeter IMF and a 9-50 $M_{\odot}$ range for progenitors of NS we have $k_{\alpha}$ = 0.0059 neutron star progenitors per $M_{\odot}$ of star-forming gas. The parameter $\alpha_{M\!N\!S}$ has been chosen in order to reproduce the current rate of GW events per unitary volume of Universe as found by Abbott et al. (2017a), $1540$ $yr^{-1}Gpc^{-3}$. The values for $\alpha_{M\!N\!S}$ are reported in table 1 for the different DTDs, alongside with the fractions of MNS systems merged after 20 Myr (third column) and 100 Myr (fourth column) after an initial single star formation burst. Together with figure 4 this allows us to make a simple comparison among the timescales of the different distributions.

In this way we have computed a redshift distribution of MNS rate directly comparable with the derived redshift distributions of SGRBs of G16 and ZW18. Both their relations have been multiplied by a suitable factor to reproduce the current GW rate, in order to make a comparison possible.

The results can be seen in figure 5. The rate of SGRB proposed by G16 is best represented by our DTD with $\beta=-1.5$, which means that a bottom heavy distribution (i.e. with many systems with small initial separations) is a good candidate. It is worth noting that such a value for $\beta$ gives rise to a distribution of time delays which scales as $\propto t^{-1.125}$. Lower values for $\beta$ imply higher and earlier maxima, so shorter time delays are preferred. However too short timescales, like those produced by the simple power law $\propto t^{-2}$ and the fixed 10 Myr total delay, do not produce a good agreement with the observations (see panel (b) of figure 5: the maxima appear to exceed by a large factor the G16 distributions.)

The rate distribution proposed by ZW18, on the other hand, cannot be matched by any of the considered distributions. It requires a very top heavy distribution for the initial separations (with a $\beta=6.0$) and a low probability of forming close binaries ($\alpha_{M\!N\!S}$ equal to $0.39\times 10^{-2}$), but such a model is in strong contrast with the local data on chemical evolution, as it will be shown in the next section.

4 Chemical evolution model for the Milky Way

A chemical evolution model tracks the abundance of a given element $i$ in the gas of a galaxy at different times. The ingredients and the general theory of such a simulation can be found e.g. in Matteucci (2012). In particular, we have adopted a scheme similar to those described in Matteucci et al. (2014) that will be summarized below.

The evolution of the surface gas density of the element $i$, $G_{i}(r,t)$, is described by the following equation:

[TABLE]

The first term is the product of the star formation rate (SFR) $\psi(r,t)$ and the elemental fraction $Z_{i}(r,t)$ in the gas and represents the gas removed by the ongoing star formation. The second term represents the fraction of gas restored as element i by stars born at $t-\tau(m)$ and dying at the time $t$. $\tau(m)$ is the stellar lifetime for a star of mass $m$. In particular, $\phi(m)$ is the initial mass function (IMF) and $Q_{mi}(t-\tau(m))$ is the production matrix for a star of mass $m$ as defined by Talbot & Arnett (1973). $M_{U}$ is the upper limit for the mass of a star, that we have chosen to be 100 $M_{\odot}$, while $m(t)$ is the minimum mass for a star dying at the time $t$. The mass dying at the present time is 0.8 $M_{\odot}$. The third term represents the infall of gas from outside the Galaxy. This infall follows a given rate in time at each galactocentric radius $A(r,t)$ and has a primordial (i.e. given by the Big Bang nucleosynthesis) composition, which for element $i$ is $Z_{i,0}$.

Our model follows the evolution of 31 chemical species, from H to Eu. The Milky Way is divided in concentric rings 2 kpc wide and eq. (19) is integrated on a variable time step, whose minimum is 2 Myr. The output quantities represents the averages (for chemical abundances) or the totals (for rates) in the simulated Galaxy. Rates are calculated for MNS and SNeIa alike as in eq. (1). We adopted the SFR by Kennicutt (1998), the IMF of Kroupa et al. (1993) and a double-exponential infall law (the so-called two infall model of Chiappini et al., 1997) aimed at reproducing the present surface mass densities in the halo and the disk as observed by Kuijken & Gilmore (1991), with timescales from Romano et al. (2000). No outflow is included and the instantaneous mixing approximation is retained (because their delaying effect on the chemical evolution seems negligible, see e.g. Spitoni et al., 2009). The yields are taken from Karakas (2010) for low and intermediate mass stars, from Doherty et al. (2014a, b) for super-AGB stars, from Nomoto et al. (2013) for SNeII, from Iwamoto et al. (1999) for SNeIa and from José & Hernanz (1998) for nova systems. The stellar lifetimes are those of Schaller et al. (1992).

5 Comparison with the Galactic abundances of europium and iron

In our calculation we need to specify the following parameters:

1. [i.]

the DTD of the MNS; 3. 2.

the DTD of the SNeIa; 4. 3.

the production of Eu by CC-SNe; 5. 4.

the fraction of neutron stars which produce a MNS, introduced in section 2 ($\alpha_{M\!N\!S}$); 6. 5.

the Eu produced per merging event.

We have varied (i) and (ii) according to the options proposed in section 2. Concerning (iii), we have run two sets of scenarios, one in which CC-SNe are allowed to produce Eu and the other in which they are not. Finally, (iv) and (v) have been chosen according to the following observational constraints:

• •

$\alpha_{M\!N\!S}$ has been fixed to the value found in section 3 to reproduce the cosmic LIGO/Virgo MNS rate (Abbott et al., 2017a);

• •

the yield of Eu has been fixed to reproduce the solar absolute abundance of Eu observed by Lodders et al. (2009), in particular comparing it with the simulated abundance of Eu present in the Galaxy 9 Gyr after the Big Bang (which is about the age of formation of the Solar System). Given that different values of $\alpha_{M\!N\!S}$ produce different predicted Eu abundances in the Sun, we have chosen a value that reproduces the observations in the case of intermediate timescales (models 4A and 4B in table 2) and we have kept it constant. This choice makes easy to see the effects of different timescales on the total amount of Eu in the solar neighborhood.

5.1 Eu production sites

As discussed in the Introduction, there are two possible sites for Eu production: MNS and CC-SNe. In this paper, we have tested two different scenarios: one in which Eu is entirely produced by MNS and another in which Eu is produced both by MNS and CC-SNe.

MNS, as a production site of Eu, have been explored thoroughly by several authors (e.g. Korobkin et al., 2012) that proposed a yield in the range $10^{-7}-10^{-5}$ $M_{\odot}$ of Eu per event. More recently, the observation of the kilonova AT2017gfo has allowed us to estimate the yield of Eu in the range $(3-15)\times 10^{-6}$ $M_{\odot}$ of Eu per event (Evans et al., 2017; Tanvir et al., 2017; Troja et al., 2017). We have chosen to represent MNS as systems of two 1.4 $M_{\odot}$ neutron stars with progenitors in the 9-50 $M_{\odot}$ mass range.

We have chosen an empirical yield of Eu equal to $4.0\times 10^{-6}$ $M_{\odot}$ when the MNS are the sole producers, and $1.5\times 10^{-6}$ $M_{\odot}$ when Eu is co-produced by CC-SNe. This second value is slightly lower than the estimate from the kilonova AT2017gfo but well inside the theoretical range.

CC-SNe, on the other hand, have a less clear role in the r-process elements production. In this paper, we have used a slightly modified version of the provisions found in Argast et al. (2004), model SN2050, and used also in the paper Matteucci et al. (2014):

• •

a constant yield of $3.8\times 10^{-8}$ $M_{\odot}$  of Eu for 20-23 $M_{\odot}$  mass range stars;

• •

a decreasing yield from $3.8\times 10^{-8}$ $M_{\odot}$  of Eu for a 23 $M_{\odot}$  star to $1.7\times 10^{-9}$ $M_{\odot}$  of Eu for a 50 $M_{\odot}$  star.

We do not consider any dependence on metallicity of these yields, nor do we take into account the distribution in rotational velocities of high-mass stars and its effect on the yield of Eu.

If the timescale of Eu production were short, we would expect to see a plateau at low metallicities ($\text{[Fe/H]}<-1.0$) followed by a decline as the Fe abundance grows. This trend is found in the observational data (see figure 6), and can be explained supposing that the main producers of Fe, namely the SNeIa, have lower mass progenitors than CC-SNe and longer gravitational delay times than MNS. Therefore, the bulk of Fe is released to the ISM on a longer timescale with respect to Eu. Initially, Eu and Fe are produced at a similar rate by only massive stars and this generates the plateau. Then, SNeIa start to explode, producing Fe but not Eu, progressively decreasing the [Eu/Fe]. This behaviour is similar to that of $\alpha$-elements (like O, Mg, Ca and Si) that are likewise produced in high mass short-lived stars, and is known as time delay model (see,e.g. Matteucci, 2012).

5.2 Conversion of the cosmic rate to a Galactic rate

As a consistency check, we have developed a simple conversion procedure to infer the Galactic MNS rate from from the cosmic LIGO/Virgo rate of Abbott et al. (2017a). We have considered the total luminosity of a unitary volume of the Universe as derived by the integration of the Press-Schechter function, equal to $2\times 10^{17}\times h$ $L_{\odot}$, as reported by Mo et al. (2010). Here, $h$ is the dimensionless parametrization of the Hubble constant in units of 100 $\text{km}\,\text{s}^{-1}\,\text{Mpc}^{-1}$.

Considering a total baryonic mass for our Galaxy equal to $(5-10)\times 10^{10}$ (as found by e.g. Sofue, 2013; McMillan, 2016) and a mass-to-luminosity ratio for a spiral galaxy equal to 5 (Mo et al., 2010), we have estimated the total luminosity of the Milky Way. The result is that there are $\sim(1-2)\times 10^{7}\times h$ Milky Way-equivalent (MWe) galaxies per $\text{Gpc}^{3}$. Choosing a value for $h$ equal to 0.7 as most cosmological observations suggest (Planck Collaboration et al., 2016) produces a central value for the obtained interval of $\sim 1\times 10^{7}$ $\text{MWe}\,\,\text{Gpc}^{-3}$. Therefore, we expect that the galactic rate of MNS is $\simeq 10^{-7}$ times lower than the Abbott et al. (2017a) value, or 1.54_{-1.22}^{+3.2}\times 10^{-4}$$\text{yr}^{-1}. This figure compares well to previous Galactic estimations based on binary pulsars such as $0.83_{-0.7}^{+2.1}\times 10^{-4}$ $\text{yr}^{-1}$ (Kalogera et al., 2004). In our chemical evolution model the current rate of MNS is not a free parameter, but results from the convolution of the star formation history of the Milky Way with the DTD, modulo the fraction $\alpha_{NMS}$ derived from the cosmic constraint. We will compare the model results to these figures.

5.3 Simulations and results

We ran several models, as reported in table 2. In the first column we put the name of the model, in the second column is indicated if the CC-SNe have been considered as production site for Eu or not and in the third and fourth columns are indicated the DTDs used, respectively, for SNeIa and MNS. In the fifth and sixth columns are reported the values for the occurrence probabilities of SNeIa and MNS. The occurrence probabilities $\alpha_{Ia}$ have been chose so as to obtain the estimated current rate of Type Ia supernovae, $1.8\times 10^{-3}\,S\!N\,yr^{-1}$ in the Milky Way (Li et al., 2011), while $\alpha_{MNS}$ has been fixed by the fit of the local volumetric SGRB rate (see Sect. 3). In the seventh column is reported the adopted Eu yield per MNS event, in the eighth column the predicted current rate of MNS in the Milky Way and in the ninth and tenth columns are reported the predicted absolute solar abundances of Fe and Eu.

As it is possible to see in table 2, the current galactic MNS rate turns out smaller than the estimate in Sect. 5.2 obtained extrapolating the Abbott et al. (2017a) value, but it appears to be in very good agreement with Kalogera et al. (2004) estimate. Given that Galactic and cosmic star formation histories are different we do expect some differences in the current MNS rate, at fixed DTD and $\alpha_{MNS}$. The effect of the Galactic SFR on the MNS rate can be clearly seen in figure 7, where the computed MNS rates are shown: at around 1 Gyr from the formation of the Galaxy, when the first infall episode terminates, the merging rate rapidly decreases. The magnitude of such a decrease depends on the timescale of the adopted DTD: longer ones act as smoothing masks. The rapid oscillations at very early times are caused by the threshold in the SFR adopted in the chemical model, as described in Sect. 4.

As a first step, it has been verified the ability of different DTDs to reproduce the solar abundances of Eu and Fe.

For the Sun, Lodders et al. (2009) have determined:

[TABLE]

[TABLE]

These values should be compared with the predicted ones in table 2.

We have kept the same yields of Fe for SNeII and SNeIa during all the tests and instead varied the DTD, using the two possible distributions listed in section 5.1. The values reported in Table 2 show that the solar Fe abundance in our model is in excellent agreement with the observational value. For the Eu contribution from MNS, the yields were actually tuned to reproduce the solar abundance with our model. We can see that the required production is quite compatible with theoretical estimates, as well as the empirical value from of the kilonova, as mentioned in Sect. 5.1.

Now we can turn to study the [Eu/Fe] vs [Fe/H] relation. We remind that the [X/Y] notation refers to the logarithm of the ratio between the two elements X and Y, normalized with respect to the same logarithmic ratio in the Sun. In figure 8 we show the effects of different DTDs on the [Eu/Fe] vs [Fe/H] relation when CC-SNe are not allowed to produce Eu. In particular, we are interested in reproducing the decreasing trend of [Eu/Fe] vs [Fe/H] in the disk stars (i.e. those at $\text{[Fe/H]}>-1.0$). We will not study the spread in the datapoints in the halo stars as done by other authors (e.g. Cescutti et al., 2015; Wehmeyer et al., 2015), since our model is homogeneous. This spread is likely due to local polluting events.

We can categorize the models in three groups: those that totally fail to produce the decreasing trend in the disk stars, those that produce an insufficient decreasing trend and those that produce an acceptable decreasing trend. The models in the first category usually fails also to reproduce the low metallicity plateau. As it is possible to see in panels (a) and (c), the DTD derived in this paper falls in the first category, meaning that coalescence timescales longer than 300 Myr are unable to reproduce the [Eu/Fe] vs [Fe/H] relation along the entire range of metallicities.

On the other hand, three out of four models shown in panels (b) and (d) reproduce the decreasing trend at high metallicity. The models that produce the best correspondence are the 7A/B and the 8A/B. In particular, the 7A/B also reproduce the plateau at $\text{[Fe/H]}<-1.0$ and can be selected as the best model. Therefore, we can say that a DTD $\propto t^{-2}$ is adequate to reproduce the [Eu/Fe] abundance pattern in the Milky Way stars, if MNS are the only contributors to Eu.

In figure 9 we show the models where CC-SNe co-produce Eu alongside MNS. The impact of the CC-SNe contribution is very large since with the prescriptions used here no less than 60% of the solar Eu should come from CC-SNe. Here, all the DTDs produce a decreasing trend in the [Eu/Fe] vs [Fe/H] relation. In particular, DTDs with longer timescales (like the one derived in this paper) become indistinguishable from one another (see panel a). For shorter timescales (panel b) there are slight differences between the models. One more time, shorter DTDs (7BS, 8BS) are favoured over longer ones (5BS). The low metallicity plateau, although deformed by the knee at $\text{[Fe/H]}\sim-3.3$, is visible nonetheless. Such a knee is caused by the death of massive stars in the lower region ([20-23] $M_{\odot}$) of the allowed mass range, where the yields of Eu are higher. We do not show the graphs from the 1AS to the 8AS model because they are substantially equal to those shown in figure 9.

In figure 10 we show in detail the effect of changing DTD for SNeIa (panel a) and activating/deactivating Eu production in CC-SNe (panel b) in the abundance range [Fe/H] > -1.5. Using the G05 DTD for SNeIa (panel a, dashed lines), that produces longer timescales than $\propto t^{-1}$, slightly increase the steepness in $\text{[Fe/H]}<-0.5$, while has substantially no impact in the $\text{[Fe/H]}>-0.5$. However, it is not sufficient to make our DTD (with $\beta=-1.5$) acceptable. Instead, the DTD $\propto t^{-2}$ offers an optimal profile in both cases (models 7A and 7B). If we observe panel (b), we can see that the activation of Eu production in CC-SNe strongly improves the prediction when using our DTD (model 4BS), even if it remains worse than in the case of shorter timescales.

Finally, in figure 11 we show the plots relative to two popular choices for MNS: the DTD $\propto t^{-1.5}$ and the DTD $\propto t^{-1}$. We recall from section 4 that, even if they are not the preferred models to reproduce the SGRB redshift distribution, they cannot be entirely ruled out. As we can see, we confirm what has been found by other studies (e.g. Côté et al., 2017a): a DTD $\propto t^{-1}$ is unable to recover the decreasing trend in the thin disk abundance pattern (panel a). If we activate the Eu production in CC-SNe (panel b), we obtain a relation that, although decreasing, is not sufficiently steep. A more satisfactory prediction is given by the DTD $\propto t^{-1.5}$ (models 6A/B and 6AS/BS), both with and without CC-SNe producing Eu. However, the steepness of the curve is one more time smaller than the observed one.

5.4 A variable $\alpha_{M\!N\!S}$

CC-SNe seem to provide a way to (partially) reconcile local and cosmic data. Another possible choice is to relax the assumption of constancy of the occurrence probability of close binary neutron stars $\alpha_{M\!N\!S}$. This parameter depends on the physics of the stellar formation and the efficiency of the common envelope phase in the late stages of stellar life. Population synthesis models (e.g Giacobbo & Mapelli, 2018b) indicate that metallicity plays an important role. Therefore, we tested a variation of the MNS occurrence frequency in the best fit model for the SGRB redshift distribution (the 4A, i.e. the one that uses the DTD derived in this paper with $\beta=-1.5$). In particular, we assumed that $\alpha_{M\!N\!S}$ depends on [Fe/H], called $Z_{F\!e}$ from now on. All the other parameters are left the same, as reported in table 2, and the value for $\alpha_{M\!N\!S}$ reported in the table for the model 4A will be indicated as $\tilde{\alpha}_{M\!N\!S}$. We tested two dependencies of $\alpha_{M\!N\!S}$ on $Z_{F\!e}$:

1. [i.]

model 4AV1, where $\alpha_{M\!N\!S}$ varies as:

[TABLE] 3. 2.

model 4AV2, where $\alpha_{M\!N\!S}$ varies as:

[TABLE]

The evolution of the parameter $\alpha_{M\!N\!S}$ is reported in figure 12. For the primordial gas, the absolute abundance of iron has been set to zero, and the initial value for $\alpha_{M\!N\!S}$ of 27.5 % has been chosen in order to best fit the observational data. Then, the sharp decrease of $\alpha_{M\!N\!S}$ in the model 4AV2 (black line) is produced by the rapid increase in metallicity caused by the death of the very first stars. This procedure is somewhat arbitrary but it is useful to show which variation of $\alpha_{M\!N\!S}$ can reconcile the SGRB cosmic rate with the [Eu/Fe] vs. [Fe/H].

We have verified the ability to reproduce the current GW rate in the Milky Way and the absolute Eu abundance in the Sun. It happens that for both models 4AV1 and 4AV2, the predicted current MNS rate differs less than 0.1% from the value reported in table 2. Speaking about the absolute solar Eu abundance, the model 4AV1 predicts a value equal to $3.79\times 10^{-10}$ and the model 4AV2 predicts a value equal to $3.88\times 10^{-10}$. Therefore, we can conclude that the impact of this enhanced early occurrence on the current MNS rate and on the total amount of Eu in the Galaxy is negligible. On the other hand, a variable $\alpha_{M\!N\!S}$ deeply influences the [Eu/Fe] vs [Fe/H] relation, as it is possible to see in figure 14. In fact, we are able to recover the decreasing trend in the disk stars, although the shape of the curve is not optimal (panel b). The model 4AV2 is also able to reproduce the low-metallicity plateau.

Finally, we have verified the impact of this new hypothesis on the cosmic rate of MNS. In order to do this, we have supposed two evolutionary patterns for the metallicity at high redshift. In the first one, the [Fe/H] evolves as predicted in the chemical evolution models of elliptical galaxies (De Masi et al., 2018). We have supposed a Salpeter (1955) IMF, an exponential infall on timescales of 0.5 Gyr and a SFR efficiency of 20. These parameters produce a transition between primordial [Fe/H] values to [Fe/H]$\sim$-1.0 dex in 80 Myr. Instead in the second one, the [Fe/H] evolves as predicted for the Milky Way, taken as representative of the spiral galaxies, on a timescale four times larger. The real evolution of the cosmic metallicity is expected to be intermediate between these two extreme cases. The results are shown in figure 13. The metallicity-dependent models differ from the non-metallicity dependent one only at very high redshift. Moreover, at least in case of a rapid chemical evolution, the predicted MNS rate does not differ too much neither at very early times. This fact demonstrates that a metallicity dependent $\alpha_{M\!N\!S}$ does not invalidate our derived DTD with respect to the cosmic MNS rate.

6 Conclusions

In this paper, we have derived a new DTD for MNS starting from the initial separation and mass distributions of neutron star binary systems. We have assumed a flat distribution of NS masses and a power law distribution for the initial separations, with exponent $\beta$. We have tested four different values for $\beta$: 0.9, 0.0, -0.9 and -1.5, similar to what suggested for the DTD of SNeIa (Matteucci et al., 2009) as derived from the double degenerate (DD) model of Greggio (2005).

We have tested our assumptions on the DTD of MNS and SNeIa on the following constraints:

1. [i.]

the cosmic SGRB rate; 3. 2.

the evolution of Eu abundance in the Galaxy, in particular the [Eu/Fe] vs [Fe/H] relation.

In this context, we have also tested the effect of core-collapse SNe as producers of Eu, besides MNS. Finally, we have tested a metallicity-dependent occurrence probability for MNS (the parameter $\alpha_{M\!N\!S}$). Our main results can be summarized as follows:

1. [i.]

Our derived DTD for MNS with $\beta$=-1.5 provides the best fit to the cosmic SGRB rate (as found by Ghirlanda et al., 2016). Also our DTD with $\beta$=-0.9 and the DTD going like $\propto t^{-1}$ give a good fit. With these three DTDs, the average timescale of coalescence of MNS is 300-500 Myr. Shorter timescales produce too many events at high redshift (2-2.5), while longer timescales produce too few events over the entire range of redshift. This result is in agreement with the one found by Fong et al. (2017) that suggest a DTD $\propto t^{-1}$ for MNS to explain the proportion of SGRB incoming from early-type galaxies. 3. 2.

The other redshift distribution found by Zhang & Wang (2017) can be reproduced only with a very top heavy (i.e with a lot of binaries with large separations) initial separation distribution, that corresponds to use our DTD with $\beta=6.0$. By the way, this extreme DTD cannot explain the Eu abundances in the Milky Way, independently from the contribution of CC-SNe to the Galactic Eu abundance. 4. 3.

When the MNS are the sole producers of Eu in the Galaxy, the models which better reproduce the [Eu/Fe] vs [Fe/H] pattern are those with short (< 30 Myr) coalescence timescales, in particular the constant total delay (stellar lifetime plus coalescence time) of 10 Myr and the $\propto t^{-2}$ distribution, in agreement with Matteucci et al. (2014) and Côté et al. (2018a). On the other hand, ourderived DTD is unable to reproduce the decreasing trend of [Eu/Fe] vs [Fe/H] observed in disk stars for reasonable values of the parameter $\beta$. 5. 4.

When the MNS are the sole producers of Eu in the Galaxy, their yield should be no more than $4.0\times 10^{-6}$ $M_{\odot}$ per event. Higher yields overestimate the absolute Eu abundance in the Sun. This is in agreement with the lower end of the estimations of Evans et al. (2017) and Troja et al. (2017) in the kilonova AT2017gfo, but much higher than the estimate of Smartt et al. (2017). 6. 5.

When CC-SNe co-produce Eu under the prescriptions of Argast et al. (2004) (their model SN2050), they become the main production site and dominate the [Eu/Fe] vs [Fe/H] relation. This can reconcile the short timescale required to explain the Eu abundance in the Milky Way with the longer timescale required to explain the SGRB redshift distribution. However, in this case, the yield of Eu per MNS event is reduced to $1.5\times 10^{-6}$ and falls (slightly) below the range of value found by Evans et al. (2017) and Troja et al. (2017). 7. 6.

The influence of different DTDs for SNeIa (at least of those tested here) on [Eu/Fe] vs [Fe/H] relation seems mostly negligible, and it is of the order of $\sim 10\%$. 8. 7.

The occurrence probability $\alpha_{M\!N\!S}$ found for SGRBs produces a Galactic rate similar to that found by Kalogera et al. (2004) from pulsar luminosities (of the order of $1-2\times 10^{-2}$). It turns out $\alpha_{M\!N\!S}\sim 10^{-2}$, interestingly close to the occurrence probability of SNeIa as determined form the cosmic rate of SNIa in the local universe (Cappellaro et al., 2015). 9. 8.

A metallicity-dependent occurrence probability, $\alpha_{M\!N\!S}$, can increase the production of Eu at early times and at low [Fe/H] values ($\leq-1.0$), enabling a DTD with average coalescence timescales of 300-500 Myr to reproduce both the SGRB redshift distribution and the [Eu/Fe] vs [Fe/H] relation in the Milky Way, when the MNS are the sole Eu producers.

In general, allowing for the CC-SNe to form Eu helps reconciling the necessity of a short timescale in Eu production with the long timescales required for the SGRBs. However, this choice makes single massive stars the main Eu production site, by contributing no less than $\sim$60% to the total Eu present in the Sun, and this fact is in contrast with the current understanding about heavy r-process elements nucleosynthesis. On the other hand, a time decreasing occurrence probability of MNS provides instead the correct production rate of Eu, without the contribution of CC-SNe, but this assumption is somewhat arbitrary and needs to be checked in the future by means of detailed population synthesis models.

Acknowledgements

PS acknowledges Carlo De Masi and Benoit Côté for some useful suggestions. LG thanks Riccado Ciolfi for stimulating scientific discussions. FM acknowledges funds from University of Trieste (FRA2016). GC acknowledges financial support from the European Union Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 664931. This work has been partially supported by the the EU COST Action CA16117 (ChETEC).