The Gravitational waves merger time distribution of binary neutron star systems
Paz Beniamini Tsvi Piran

	
	Merger time of binary neutron stars]The Gravitational waves merger time distribution of binary neutron star systems
	
P. Beniamini & T. Piran]
Paz Beniamini^1,2E-mail: paz.beniamini@gmail.com
& Tsvi Piran^3



^1Department of Physics, The George Washington University, Washington, DC 20052, USA 

^2Astronomy, Physics and Statistics Institute of Sciences (APSIS)
	
^3Racah Institute of Physics, the Hebrew University, 91904, Jerusalem, Israel
	

	

	
[
    [
    December 30, 2023
=====================


	Binary neutron stars (BNS) mergers are prime sites for r-process nucleosynthesis. Their rate determines the chemical evolution of heavy elements in the Milky Way. 
	The merger rate of BNS is a convolution of their birth rate and the gravitational radiation spiral-in delay time. Using the observed population of Galactic BNS we show here that the lifetimes of pulsars in observed BNSs are sufficiently short that the ages of BNSs have little to no effect on the observed merger time distribution. We find that at late times (t≳ 1 Gyr) the gravitational wave delay time distribution (DTD) follows the expected t^-1. However, a significant excess of rapidly merging systems (between 40-60% of the entire population) is apparent at shorter times.
	Although the exact shape of the DTD cannot be  determined with the existing data, in all models that adequately describe the data we find at least 40% of BNSs with merger times less than 1Gyr.
	This population of rapid mergers implies a declining deposition rate of r-process materials that is consistent with several independent observations of heavy element abundances in the Milky Way.
	At the same time this population that requires initial binary separations of roughly one solar radius clearly indicates that these binaries had common envelope progenitors. Our results suggest that a significant fraction of future LIGO/Virgo BNS mergers would reside in star forming galaxies.  



	gravitational waves – stars: neutron – Galaxy: abundances – stars: evolution 




§ INTRODUCTION


The merger of binary neutron star (BNS) systems results in a blast of gravitational wave (GW) radiation, in a short gamma-ray burst (GRB) and in the synthesis of r-process elements at a rate that can account for the overall formation of those elements in the Universe. The connection between these different aspects of BNS mergers  <cit.> has recently been given very strong support due to the detection of GWs from GW170817 which were accompanied  by a GRB <cit.>, a macronova (powered by the radioactive decay of r-process elements) and a multiwavelength afterglow. 
Interestingly the implied rate and yield of heavy r-process element is consistent with that required for the production of heavy r-process elements in the Galaxy <cit.>.


The rate of BNS mergers depends on the birth rate of these systems as well as on the delay times from birth to merger, commonly known as the `delay time distribution' (DTD). As neutron stars' progenitors are massive stars their birth rate follows the star formation rate (SFR) with a minimal delay. The DTD is dominated by this GW spiral-in time. As a BNS loses energy due to GW radiation its orbit shrinks and eventually the neutron stars merge. The time until merger t_ m∝ a^4 (1-e^2)^7/2 depends on the initial semi major axis, a, and the eccentricity, e, of the BNS. As such the determination of the DTD and in particular of the minimal time delay, t_ min will provide valuable information on the unknown final stages of evolution of BNS progenitors. 
Here we explore the DTD  using the observed Galactic population of BNS systems. This approach is complementary to attempts to explore the DTD by comparing the rate of short GRBs (believed to be associated with BNS mergers) and the SFR <cit.>. 
An independent single observation arises from GW170817. The host galaxy of this event is an S0-type galaxy with a very low star formation rate, suggesting a time delay between the BNS formation and merger of ≈ 1-10 Gyr.


The GW merger time, t_ m, depends strongly on a, and since many BNSs exhibit low orbital eccentricities, it is insightful to consider the expected DTD resulting from an initial power-law distribution of semi major axes, dN/da∝ a^-n with circular orbits (e=0). The result is dN/dt_ m∝ t_ m^-(n+3)/4. In particular, if the semi-major axis is uniformly distributed in log-space, dN/dlog a=const, one gets the well known result dN/dt_ m∝ t_ m^-1 <cit.>. Note also that due to the strong dependence of t_ m on a, the distribution of semi-major axes must be very steep in order to get significantly steeper dN/dt_ m profiles. 


We use the observed distribution of merger times of the Galactic BNSs to derive the DTD in a model-independent way, directly from observations. The paper is structured as follows. In <ref> we describe the observed sample. We then study analytically the connection between the distribution of delay times at birth and the merger time (and age) distribution viewed at some later time in <ref>. In <ref> we perform a likelihood analysis in order to determine the underlying DTD from the observed population. We examine some implications on the total number of un-merged systems and observable pulsars remaining in BNSs in the Galaxy today in <ref>. We compare our results with other findings in <ref> and we summarize and discuss some broader implications in <ref>.



§ GALACTIC BNS OBSERVATIONS


Our sample includes the 15 observed field Galactic BNS systems. We exclude two Galactic BNSs that reside in globular clusters. Table <ref> provides the semi-major axes, eccentricities, merger times and spin-down times of these systems. The observed merger time distribution is shown in figure <ref>. 

Clearly this  sample  does not include all Galactic BNSs, as the observed systems are only those that have an active pulsar pointing towards us. Since pulsars live a finite time there are many systems that cannot be observed. The BNS pulsars are mostly recycled ones <cit.> and their true ages[We consider as the true age of the system the time since the formation of the second neutron star.] may be considerably younger than  their characteristic dipole spin-down times, τ=P/2Ṗ, <cit.>. 
This is because for any pulsar with a braking index n>1, the ratio τ=P/(n-1)Ṗ∝ P^n-1 increases with P, and always remain an upper limit on the true age <cit.>.
For example, in the double pulsar system, J0737-3039, the two pulsars (which are expected to be of a similar age) have spin-down ages of 210 and 50 Myr <cit.>. This is much longer than what is expected to be their true age difference. 

 



§ GENERAL CONSIDERATIONS




 §.§ The delay time distribution - DTD


The merger time distribution, dN/dt, is a convolution of the intrinsic DTD, D(t), with the BNS birth rate, r(t), where time here is measured backwards from the present day (i.e. today we have t=0 and at earlier times, t is larger): 

    dN/dt=∫_0^t_ H r(t')D(t'+t)dt',

where t_ H is the Hubble time. We will consider DTDs that have a minimal t_ min and maximal t_ max range. 

The distribution dN/dt is different from the observed merger time distribution  dN_ obs/dt since, as mentioned earlier, the observed distribution depends on detection of the corresponding pulsars and those have a finite lifetime. 
To account for the finite lifetimes of the pulsars in the BNS systems we introduce a typical lifetime t_ a in which one of the pulsars in the BNS system can still be observed.  In reality, there is a distribution of lifetimes. As there is not enough information today to reveal this distribution we will take t_ a to  be a constant. If there is a distribution of lifetimes than t_ a,min < t_ a < t_ a,max the transition between the regimes discussed below will be gradual, but the basic qualitative features will remain. The regime t≪ t_ a (t≫ t_ a) should then be interpreted as t≪ t_ a,min (t≫ t_ a,max). We show later (see <ref>) that the value of t_ a is likely small compared to the typical delay times. However, for the sake of generality, we make no assumptions on the value of t_ a in the following discussion.



We write the number of observed BNSs with a given time to merger, dN_ obs/dt, as:

    dN_ obs/dt=∫_0^t_ ar(t')D(t'+t)dt'.

Since, as we show later, the life times of the observed BNSs as active pulsars is much shorter than the scale on which the Galactic star formation rate (SFR) changes  we use here a constant Galactic SFR, r(t)=r_0. This simple assumption has been suggested in the literature as a good approximation for the star formation rate of the Milky Way (MW) during the majority of its lifetime <cit.>.  For r(t)=r_0  and for t_ a<t_ min, the observed merger distribution is:

    dN_ obs/dt∝{[                     0         t<t_ min-t_ a;                 const t_ min-t_ a<t≪ t_ min;                  D(t)            t≫ t_ min. ].

Note that for t_ a≪ t_ min < t the observed distribution is simply proportional to the DTD:  dN_ obs/dt∝ D(t) regardless of the functional form of D(t). Consider now the opposite limit, t_ a>t_ min (assuming once more a constant formation rate). For t> t_ a, that is for 
merger times much longer than the age of the systems, the pulsars stop shining on a time scale much shorter than their merger time. Hence, dN_ obs/dt ∝ D(t) as before. 
But, of course, there are many more (t_ H/t_ a) BNS systems that we cannot observe. Alternatively, at much shorter times, t≪ t_ a the observed population is reduced compared to an extrapolation from the t≫ t_ a regime since only a fraction t/t_ a of the systems born in the last t_ a years with a delay time t have not merged yet. dN_ obs/dt is therefore proportional to (t/t_ a)D(t). At delay times even lower than t_ min, the only systems observable are the tail of the systems born with a delay of ∼ t_ min (which by construction are the type of systems that dominate the intrinsic population) that were formed t_ min-t years ago. For t≪ t_ min the two latter factors are constant, and thus dN_ obs/dt is constant.

Consider, for example, a power-law DTD:

    D(t)=1-b/t_ max^1-b-t_ min^1-b  t^-b   for   t_ min≤ t ≤ t_ max .

For  an index b> 1 and t_ a>t_ min we obtain

    dN_ obs/dt∝{[            const        t≪ t_ min;            t^1-b  t_ min≪ t≪ t_ a;             t^-b t_ a≪ t≪ t_ max. ].

Overall we find that the finite lifetime of pulsars implies that for merger durations longer than t_ a the observed merger time distribution simply follow the DTD. For durations shorter than t_ a we observe fewer systems than what the DTD would predict.



 §.§ The Galactic BNS distribution
 
The observed distribution of Galactic BNS merger times is shown in figure <ref>.
At times much larger than the typical ages of the pulsar, dN_ obs/dt ∝ D(t), i.e. it is directly related to the DTD.
It is evident from figure <ref> that at t≳ 1 Gyr, dN/dlog t_ m is approximately constant, suggesting D(t) ∝ t^-1. However, at lower merger times we find an observed excess of systems with 0.05 ≲ t_ m≲ 1 Gyr (seven out of fifteen systems, representing a fraction Z_<1 Gyr= 0.47 of the population), as compared with the D(t)∝ t^-1 model predictions (for which  Z_<1 Gyr=0.19, as we discuss later in <ref>). At the same time, as we argued above, the finite lifetime of the pulsars leads instead to a suppression of the number of systems with t_ m < t_ a. 


D(t)∝ t^-1 is a good description of the observations for t>1 Gyr. The excess of short duration mergers (relative to predictions of the D(t)∝ t^-1 model) at t_ m≲ 1 Gyr must be an intrinsic feature of the DTD. 
Note that the possibility  that this excess is due to an increase of the BNS formation rate in the last 1 Gyr is highly unlikely. This time-scale  is much smaller than the scale on which significant Galactic changes took place. 





 §.§ Ages of observed BNS systems


As mentioned in <ref>, the ages of BNS systems are less or equal than their spin-down times. The distribution of observed ages of BNS systems, is naturally affected by the DTD (at older ages, a larger fraction of systems have already merged). The number of systems with age t_ age is given by

    dN/dt_ age   =     r(t_ age)∫_t_ age^∞ D(t) dt 
       ∝   {[                     const             t_ age<t_ min;                t_ age^1-b t_ min < t_ age < t_ max. ].

where  we have assumed in the second line a constant formation rate and a power-law DTD with b>1.

One of the binary systems, J1906+0746, has a very short spin down time, P/2Ṗ=10^-4 Gyr. 
This value is much smaller than the smallest  observed delay time (see table <ref>).
If indeed P/2Ṗ reflects  an upper limit on the real age of this pulsar then, independently of the shape of the delay time distribution,  
equation <ref> implies that there should be many older systems.
As an illustration, there should be roughly 100 systems with P/2Ṗ<0.01 Gyr for each system with P/2Ṗ<10^-4 Gyr. Even taking into account a possible observational bias according to which older system may be harder to detect this ratio is inconsistent with the observations. 

To explore the implication of this inconsistency we consider the  P/2Ṗ values in table <ref> as upper limits on the true ages of the systems. We assume that the systems live for a fixed time t_ a  and search using a Monte Carlo simulation for a value of t_ a for which a simulated set of N ages (where N is the number of observed BNS systems) resulting from a given b,  t_min and t_a is consistent with these upper limits. Namely we look for the maximal value of t_a for which more than 5% of the trials are consistent. Larger values of t_a are ruled out at 2 σ.

For 0.01 ≤ t_ min≤ 0.05 and 1≤ b ≤ 1.5, as suggested by the observed merger times [The results hold as long as there is no significant drop in the age distribution given by Eq. <ref> above t_ min.], we find t_ a<0.03 Gyr at a 3σ confidence level. This result, which depends weakly on the choice of t_ min and b depends very strongly on the very short P/2Ṗ of  J1906+0746. 
However, even if the spin-down age of the pulsar in J1906+0746 is disregarded as an anomaly, the remaining distribution of spin-down ages still constrains t_ a<0.3 at a 2σ level. That is even if t_ a>t_ min,  the dynamical range between the two cannot be large. Therefore, the effect of t_ a on the observed distribution must be relatively small. This conclusion is  consistent with the estimate provided in <ref> based on the number of observed Galactic BNS systems in comparison with estimates of the BNS formation rate.
As discussed in <ref>, this means that the observed distribution of merger times roughly traces the observed ones, and that the exact value of t_ a is not required to deduce the shape of the DTD.

Following these results we assume in the following that t_ a<t_ m. This means that the orbits of the majority of systems don't evolve appreciably during the time  that their pulsars remain observable and hence the observed merger time distribution reflects the intrinsic one. We also verify later (see <ref>), using a maximum likelihood analysis, this assumption. 



§ LIKELIHOOD ANALYSIS


We consider several models that can be fitted to the observed merger time distribution. Our models involve one or two parameters (denoted here as X_i for the more general case) that describe the DTD.
We search for the values of X_i that maximize the likelihood of obtaining the observed distribution of Galactic BNS merger times.
To do so we construct a likelihood function as follows

    L(X_i)=Π_j  P(t_ m,j|X_i)

where j goes over the observed systems in the sample, P(t_ m,j|X_i) is the probability that a system was born with a delay time t_ m,j given the delay time distribution specified by X_i. 



 §.§ A t^-1  delay time distribution


Our first choice for the DTD is a simple power-law t^-1. This is motivated, 
as mentioned in <ref>, by the strong dependence of t_ m on the initial BNS separation  <cit.>. Moreover, an inspection by eye of the long merger time part of the distribution (which is not influenced by pulsar's lifetimes and the choice of the edges of the distribution) suggests that this model is consistent at least for long merger durations.
The observed delay times depend only weakly  on t_ max which we fix here to be t_ max=10^6 Gyr. 
The DTD is in this case a function of a single parameter t_ min. The likelihood  is maximized for t_ min=0.035 Gyr. 
Using a Kolmogorov-Smirnov (KS) test, we compare the observed merger time distributions with the model. The two distributions are shown in figure <ref>. The KS test rules out the t^-1 model at a 2σ confidence level. As mentioned in <ref>, there is a much larger fraction of short mergers in observations (Z_<1  Gyr=0.47) than expected for the best D(t)∝ t^-1 (Z_<1  Gyr=0.19). This is the main source of discrepancy between this model and observations. It is important to stress that, as we have seen in <ref>, introducing longer lifetimes will only exacerbate the problem as this leads to a decrease in the number of observed systems with short merger times (see figure <ref>).
Given this discrepancy we explore in the following several DTD models that involve a second free parameter.




 §.§ A general power-law delay time distribution


The observed merger times clearly span many orders of magnitude (see table <ref>). The simplest  functional form spanning a wide range is  a power-law.
A natural generalization is to use another power-law t^-b and let b≠ 1 be a free parameter in addition to t_ min. Due to the vast range between t_ min and t_ max (at least seven orders of magnitude; see table <ref>) even a small change in b can have a significant effect on the observed distribution. Indeed taking b=1.1, t_ min=0.035 (and  t_ a<t_ min) results in a distribution (denoted in figure <ref> as D_3(t)) that is consistent with  the current data. The main importance of this possibility is the simplicity of implementing it in analytic expressions such as in <ref>, <ref>.



 §.§ A short delay time population


At long merger times ≳1 Gyr, the DTD is well described by D(t)∝ t^-1. Motivated by the excess of events with short merger times,  we add a short merger time log-normal component to the t^-1 power-law: 

    D_1(t)=At^-1 + B   Lognormal(μ,1) t_ min<t<t_ max .

The second term is a log-normal distribution with σ=1 (implying a width that is of the same order of magnitude as the median). Based on the observed sample and on the modeling for the t^-1 and general power-law DTDs we fix here t_ min=0.035 Gyr, t_ max=10^6 Gyr. This model involves two free parameters: μ, f, where f is the ratio of events resulting from the log-normal distribution  out of the total number of events (and is proportional to B/(A+B)). 

The likelihood function for D_1(t) is maximized at exp(μ)=0.3 Gyr and f=0.4. It is depicted in figure <ref>. The resulting DTD is now statistically consistent with the data (see also figure <ref>). The model is insensitive to the exact value of t_ min. Values in the range 0.001≤ t_ min≤ 0.1 are all consistent with observations. 

Given the  small number of observed systems, the functional form of the DTD cannot be uniquely determined by the available data. As an illustration, we consider a broken power-law:

    D_2(t)=C	{[         (t/t_0)^-a    t_ min<t< t_0 ,;         (t/t_0)^-1 t_0 < t < t_ max . ].

We take, as before, t_ min=0.035 Gyr, t_ max=10^6 Gyr. The free parameters of D_2(t)  are a,t_0. This model is also consistent. The parameters  t_0=15, a=1.3, t_ min=0.035 Gyr  provide a statistically consistent fit to the data. Allowing for varying values of t_ min we find that consistent solutions can be found as long as 0.02≤ t_ min≤ 0.15. 






 §.§ Comparison between the models


To compare D_1(t),D_2(t),D_3(t) to D∝ t^-1 we perform a likelihood ratio (LR) test. The LR is defined as LR=2log(L_ max,a/L_ max,0) where L_ max,a, L_ max,0 are respectively the maximum likelihoods of some model a and of the null hypothesis (model 0). LR follows a χ-squared distribution where the number of degrees of freedom is given by the difference in free parameters between the two models (2 between D_1,D_2 and D∝ t^-1 and 1 between D_3 and D∝ t^-1). For all three models the improvement as compared with a t^-1 DTD is statistically significant (with a p-value for chance coincidence of p<0.05). 
Inspection by eye of figure <ref> suggests that the model with an extra short-time  log-normal population looks best. However, this is not statistically significant and it depends on the binning of the data. A  shared property of  all these consistent models is that they all result in a large fraction, Z_<1  Gyr≳ 0.4,  of BNS with merger times of less than 1 Gyr.

When a significantly larger number of BNS systems are observed in the future it will be possible to differentiate between more subtle details of the delay time distribution.
In particular this will enable to differentiate between the different models described above.  
A comparison of non-nested models with each other (i.e. comparing any two of D_1,D_2,D_3) is impossible with the regular LR test and usually requires a larger data-set to be  meaningful <cit.>. 
With coming detections of gravitational waves from BNS mergers, the SFR of the host galaxies will enable us to obtain independent evidence for this fast merging population. This will also enable us to distinguish between the different models described above.  



 §.§ The pulsar's lifetime 


Before concluding this section we return once more to the question of the pulsar's lifetime.  To address it we have carried out a likelihood analysis including t_ a as a free parameter. Figure <ref> depicts the likelihood projections for model D_3 but now with t_ a as a free parameter. Note that in this analysis we did not impose the conditions t_ age < P/2Ṗ, as with the very short P/2Ṗ of J1906+0746  this would have immediately imposed a very small t_ a (see <ref>). Still even without this constraint we find (see figure <ref>) that t_ a< t_ min is preferred, where the best fit value is  ≲ 0.02 Gyr and where in any case the region t_ a > 0.1 Gyr is practically excluded. 





§ TOTAL NUMBER OF SURVIVING AND MERGED GALACTIC BNS


We turn now to explore the unobserved populations: BNS systems that cannot be observed because their pulsars have died or are not pointing towards us and BNS systems that have already merged.
The fraction of merged systems depends on the formation rate history of BNSs in the MW and on the DTD:

    Ψ_ merg= 1/∫_0^t_ Hr(t) dt∫_0^t_ Hr(t) P(t_d<t)dt  ,

where P(t_d<t)=∫_0^tD(t)dt is the probability of a system born t years ago to have merged by today. A specific case of interest is that of a constant rate r(t)=∫_0^t_ Hr(t) dt/t_ H≡ r_0. 
In this case, for a power-law distribution, we have: 

    Ψ_ merg   =   t_ H^2-b+(1-b)t_ min^2-b - (2-b)t_ H t_ min^1-b/(2-b) t_ H(t_ max^1-b-t_ min^1-b)
       ≈    1-(t_ min/t_ H)^b-1≈ 0.5 ,

where in the last expression we have used the results for D_3(t) from <ref>.
We find a similar fraction of merged systems, Ψ_ merg(D_1)≈Ψ_ merg(D_2)≈ 0.55,  for the distributions, D_1(t),D_2(t).
We conclude that about half of the BNS systems formed in the history of the Galaxy have already merged.
Using a BNS formation rate within the MW of ∼ 4× 10^4^-1 <cit.>, this translates to a total of ∼ 2.7× 10^5 un-merged BNS systems in the MW.


We next turn to estimate the total number of BNSs with active pulsars that reside in the Galaxy today. As before we take the formation rate to be constant up to t_ a (see <ref>).  The total number of active systems in the MW is then estimated as

    N_ active=∫_0^t_ a r(t)P(t_d>t)dt≈
        r  t_ a{[                       1            t_ a< t_ min; (t_ a/t_ min)^1-b/(2-b)           t_ a> t_ min. ].

The observable number of systems in the MW is smaller than the number of active ones due to two factors. The first is Φ_ beam. This is the correction that arises since pulsars are detected by their beamed radio emission. For each pulsar with a period P and temporal beam width W, there are approximately P/2W actual pulsars in the Galaxy of the same type. 
Using W_50 as the beam width (where W_50 is the width of the pulsar's pulses at a 50% level), we find that for our sample Φ_ b≈⟨ P/(2W)⟩ ^-1=0.1. The other reduction is due to the incompleteness of the BNS pulsar population, Φ_ inc. For example, <cit.> have predicted that SKA could increase the number of known radio pulsars by a factor of ∼ 5-10, implying Φ_ inc<0.2. Overall, the observed current number of BNSs in the Galaxy, N_ obs≈ 20 implies N_ active=N_ obsΦ_ b^-1Φ_ inc^-1≳ 10^3. These estimates are consistent with estimates by <cit.> that calculated there should be ≈ 300 Galactic systems similar to each B1534+12 and B1913+16.
We can now estimate t_ a,

    t_ a≈N_ obs/Φ_ bΦ_ incr=0.015 (0.1/Φ_ b)(0.03/Φ_ inc) (4× 10^4/r )

This is consistent with the condition t_ a≲ t_ min suggested by the distribution of spin-down ages (see <ref>).




§ COMPARISON WITH OTHER RESULTS
 
Our main finding is an excess of rapid mergers over the expected t^-1 delay time distribution. It turns out that there are other venues that lead to this conclusion. An excess of rapid mergers has also been suggested on completely separate grounds following from observations of r-process element abundances.
Specifically, the high abundance of r-process elements in the ultra-faint dwarf (UFD) galaxy Ret-II as compared to the upper limits on those abundances in other UFDs suggests a neutron star merger origin <cit.>. At the same time, UFDs are known to have ceased their star formation within the first Gyr of their formation <cit.> requiring a rapid merger to enrich the gas before it can be reprocessed into the stars seen in those galaxies today <cit.>. More generally, if UFD-like galaxies form the main building blocks from which the MW halo population is composed, then rapid mergers are required to explain the high values of r-process abundance observed in extremely metal poor stars <cit.>. Furthermore, a significant population of rapid mergers, coupled with a declining star formation rate over the last few Gyrs implies a declining rate of mergers over the evolution of the MW. This is consistent with the observed rate of deposition of radioactive ^244Pu on Earth today which is considerably lower than its value 4.6 Gyr ago when the solar system was formed <cit.>. Finally, a population of rapid mergers was shown by <cit.> to provide an explanation for the declining rate of [Eu/Fe] as a function of [Fe/H] observed in the MW for ≳ -1 which is not trivially explained with a D∝ t^-1 DTD (see e.g. ). This result is depicted in figure <ref> in which we repeat the one zone chemical evolution model presented in 4 of <cit.> with the same model parameters (and for the case of a constant Galactic SFR) and changing only the DTD of the NS mergers. With no delays the results track closely the observed abundance patterns in Milky Way stars. However, a t^-1 DTD, clearly evolves too shallowly.

As mentioned in <ref>, using the sGRB redshift and peak flux distribution, <cit.> have found that dN/dt_ m∝ t_ m^-0.8 ± 0.25. There is some tension between this result and ours. However, <cit.> allow in their analysis for the existence of a second population that follows the SFR. They find that this population includes about 1/3 of sGRBs and interpret this as these bursts being “short duration" Collapsars <cit.>. These findings are consistent with our results if a significant fraction of these are genuine sGRBs. 

We can compare our results to the implications of the first NS merger observed in gravitational waves, GW170817. As mentioned earlier the time delay for GW170817 is estimated at t_ 170817≈ 1-10 Gyr. Thus, GW170817, can be used to put an upper limit on the rate of short mergers due to the following argument. The inferred rate of BNS mergers and the amount of r-process mass estimated from this single event are already large enough to dominate the production of r-process elements in the Universe (see e.g.  for a detailed discussion).
In fact, they are already somewhat larger than those inferred from the total r-process mass in the MW (albeit with still relatively large errors).
If the delay time distribution is significantly steeper than t^-1, then the inferred rate of overall mergers (and in turn their contribution to the total r-process mass) would have been much larger in the past (  ∝ (t_ 170817/t_ min)^b-1).  Applying a conservative limit, that the overall merger rate should not be more than 5 times larger than that inferred by GW170817, we find
b<1+log(5)/log(t_ 170817/t_ min)
This leads to b<1.5 for t_ 170817=1 Gyr and b<1.3 for t_ 170817=10 Gyr. These limits are consistent with the level of enhancement of short mergers found in the present work. 





§ DISCUSSION


We have explored the distribution of delay times (DTD) between formation and merger of binary neutron star (BNS) systems.
We have shown that the times since formation of the observed population of BNS systems, that are dictated by the lifetime of the pulsars, are sufficiently short so that the ages of the observed BNS systems have little to no effect on the observed merger time distribution. This distribution then directly traces the intrinsic DTD distribution. Using a maximal likelihood we found several distributions that fit the observed data. 


At long times (t≳ 1 Gyr), the DTD is consistent with the theoretically expected D∝ t^-1. At shorter times, however, there is a statistically significant excess of rapidly merging systems (almost half of the systems merge in less than 1 Gyr). We have shown that the finite pulsar's lifetime cannot explain this excess which must be intrinsic. We explored a few models with this characteristic.
One such model that provides a good description of the data is a DTD where  40% of the systems arise from a log-normal distribution with μ=log(0.3/), σ=1 while the other  60%  arise from a t^-1 distribution, with minimum and maximum merger times of t_ min=0.035 , t_ max=10^6. However there is not enough data to distinguish between this distribution and distributions that have other functional forms that include an excess of short merger time systems. 

The observed chemical evolution of heavy r-process elements in the Milky Way is estimated from the ratio of [Eu/Fe] as compared with  [Fe/H]. The decrease in [Eu/Fe] for ≳ -1 has been puzzling within the context of BNS as sites of r-process nucleosynthesis as it is not trivially explained with a D∝ t^-1 DTD (see e.g. ). As mentioned earlier this  population of rapid mergers    provide a natural  explanation for this puzzle  <cit.>.

The range of merger times, from t_ min≈ 0.035  (0.001-0.15)  Gyr[The values in bracket correspond to upper and lower estimates of t_ min]  to t_ max≈10^6  Gyr corresponds to initial orbital separations (for a canonical 1.4-1.4 M_⊙ binary) of   7  (3-10)· 10^10   cm to 5 · 10^12  cm. The minimal value is the most interesting as it provides an indication on the state of the system just before it collapses to form a neutron star. This distance is about  a solar radius. Since the progenitor is essentially heavier than that, it implies that the progenitor must have lost a significant fraction of its envelope just before the collapse, or that its companion was within the envelope during the collapse.

Our inferred DTD suggests that roughly half of the BNS systems ever formed  have merged by today. This could lead to ≈ 3× 10^5 unmerged BNS systems in the MW. That being said, the vast majority of those systems are not expected to host active pulsars.


The departure of the DTD from a simple t^-1 may point towards different channels of BNS formation.
Indeed <cit.> have shown that there are two classes of collapses that lead to the formation of the second NS in a BNS system. In roughly 60-70% of the systems, the collapse involves only a small amount of mass ejection (Δ M≲ 0.5 M_⊙) and a weak kick (v_k≲ 30), while in the rest of the systems the collapse involve larger ejected masses and kicks. The first type of collapse results in low eccentricity (e≈ 0.1) systems, while the second can result in much more eccentric orbits. Since there appears to be no correlation between e and t_ m (see table <ref>), it seems unlikely that the excess of rapidly merging systems is related to the nature of the kicks. An alternative possibility is that the excess is related to the initial separation at birth (see also ), which in turn may indicate a different stellar evolutionary part  involving a larger shrinkage of the orbit in a common envelope phase.

GW detectors are expected to detect many more BNS mergers in the coming years <cit.>. Although the coincident detection of GRBs are predicted to be much less common <cit.>, we may expect the detection of the isotropic Macronovae in most cases. These will provide us with further data concerning the rate and r-process yield of BNS mergers. The SFR history of the host galaxies could reveal information on the DTD. We expect that at least 50% of the hosts will show evidence of recent SFR. 
Advanced radio telescopes such as SKA would increase the sample of Galactic BNSs and this will enable us to improve the statistics explored here.



§ ACKNOWLEDGMENTS

We thank Ehud Nakar,  Bernard Schutz and Ben Shenhar for helpful discussions. TP was partially supported by an ERC advanced grant TReX and by the CHE-ISF center for excellence in Astrophysics. 
	
67
	natexlab#1#1[Abbott et al(2017)Abbott, Abbott, Abbott,
		Acernese, Ackley, Adams, Adams, Addesso, Adhikari, Adya, &
		et al.]GW170817
	Abbott B. P. et al., 2017, Physical Review Letters, 119, 161101
	
	[Argast et al(2004)Argast, Samland, Thielemann, &
		Qian]Argast2004
	Argast D., Samland M., Thielemann F.-K., Qian Y.-Z., 2004, , 416,
	997
	
	[Belczynski et al(2018)Belczynski, Bulik, Olejak,
		Chruslinska, Singh, Pol, Zdunik, O'Shaughnessy, McLaughlin,
		Lorimer, Korobkin, van den Heuvel, Davies, &
		Holz]Belczynski2018
	Belczynski K. et al., 2018, arXiv e-prints
	
	[Belczynski et al(2006)Belczynski, Perna, Bulik,
		Kalogera, Ivanova, & Lamb]Belczynski2006
	Belczynski K., Perna R., Bulik T., Kalogera V., Ivanova N., Lamb
	D. Q., 2006, , 648, 1110
	
	[Beniamini, Dvorkin & Silk(2018)Beniamini, Dvorkin, &
		Silk]Beniamini2018
	Beniamini P., Dvorkin I., Silk J., 2018, , 478, 1994
	
	[Beniamini, Hotokezaka & Piran(2016a)Beniamini,
		Hotokezaka, & Piran]BHP2016b
	Beniamini P., Hotokezaka K., Piran T., 2016a, , 829,
	L13
	
	[Beniamini, Hotokezaka & Piran(2016b)Beniamini,
		Hotokezaka, & Piran]BHP2016
	Beniamini P., Hotokezaka K., Piran T., 2016b, , 832, 149
	
	[Beniamini et al(2019)Beniamini, Hotokezaka, van der
			Horst, & Kouveliotou]Beniamini2019
	Beniamini P., Hotokezaka K., van der Horst A., Kouveliotou C., 2019,
	[Beniamini & Piran(2016)]BP2016
	Beniamini P., Piran T., 2016, , 456, 4089
	
	[Bromberg et al(2013)Bromberg, Nakar, Piran, &
		Sari]Bromberg2013
	Bromberg O., Nakar E., Piran T., Sari R., 2013, , 764, 179
	
	[Brown et al(2014)Brown, Tumlinson, Geha, Simon,
		Vargas, VandenBerg, Kirby, Kalirai, Avila, Gennaro, Ferguson,
		Muñoz, Guhathakurta, & Renzini]Brown2014
	Brown T. M. et al., 2014, , 796, 91
	
	[Burgay et al(2003)Burgay, D'Amico, Possenti,
		Manchester, Lyne, Joshi, McLaughlin, Kramer, Sarkissian,
		Camilo, Kalogera, Kim, & Lorimer]Burgay2003
	Burgay M. et al., 2003, , 426, 531
	
	[Champion et al(2005)Champion, Lorimer, McLaughlin,
		Xilouris, Arzoumanian, Freire, Lommen, Cordes, &
		Camilo]Champion2005
	Champion D. J. et al., 2005, , 363, 929
	
	[Corongiu et al(2007)Corongiu, Kramer, Stappers,
		Lyne, Jessner, Possenti, D'Amico, & Löhmer]Corongiu2007
	Corongiu A., Kramer M., Stappers B. W., Lyne A. G., Jessner A.,
	Possenti A., D'Amico N., Löhmer O., 2007, , 462, 703
	
	[Côté et al(2016)Côté, Ritter,
		O'Shea, Herwig, Pignatari, Jones, & Fryer]Cote2016
	Côté B., Ritter C., O'Shea B. W., Herwig F., Pignatari M.,
	Jones S., Fryer C. L., 2016, , 824, 82
	
	[Coward et al(2012)Coward, Howell, Piran, Stratta,
		Branchesi, Bromberg, Gendre, Burman, & Guetta]Coward2012
	Coward D. M. et al., 2012, , 425, 2668
	
	[Dietz(2011)]Dietz2011
	Dietz A., 2011, , 529, A97
	
	[Eichler et al(1989)Eichler, Livio, Piran, &
		Schramm]Eichler1989
	Eichler D., Livio M., Piran T., Schramm D. N., 1989, , 340, 126
	
	[Faucher-Giguère & Kaspi(2006)]Giguere2006
	Faucher-Giguère C.-A., Kaspi V. M., 2006, , 643, 332
	
	[Faulkner et al(2005)Faulkner, Kramer, Lyne,
		Manchester, McLaughlin, Stairs, Hobbs, Possenti, Lorimer,
		D'Amico, Camilo, & Burgay]Faulkner2005
	Faulkner A. J. et al., 2005, , 618, L119
	
	[Ferdman et al(2014)Ferdman, Stairs, Kramer,
		Janssen, Bassa, Stappers, Demorest, Cognard, Desvignes,
		Theureau, Burgay, Lyne, Manchester, & Possenti]Ferdman2014
	Ferdman R. D. et al., 2014, , 443, 2183
	
	[Guetta & Piran(2005)]Guetta2005
	Guetta D., Piran T., 2005, , 435, 421
	
	[Guetta & Piran(2006)]Guetta2006
	Guetta D., Piran T., 2006, , 453, 823
	
	[Guetta & Stella(2009)]Guetta2009
	Guetta D., Stella L., 2009, , 498, 329
	
	[Hotokezaka, Beniamini & Piran(2018)Hotokezaka, Beniamini,
		& Piran]HBP2018
	Hotokezaka K., Beniamini P., Piran T., 2018, International Journal of
	Modern Physics D, 27, 1842005
	
	[Hotokezaka, Piran & Paul(2015)Hotokezaka, Piran, &
		Paul]Hotokezaka2015
	Hotokezaka K., Piran T., Paul M., 2015, Nature Physics, 11, 1042
	
	[Janssen et al(2008)Janssen, Stappers, Kramer,
		Nice, Jessner, Cognard, & Purver]Janssen2008
	Janssen G. H., Stappers B. W., Kramer M., Nice D. J., Jessner A.,
	Cognard I., Purver M. B., 2008, , 490, 753
	
	[Kalogera, Kim & Lorimer(2003)Kalogera, Kim, &
		Lorimer]Kalogera2003
	Kalogera V., Kim C., Lorimer D. R., 2003, in Astronomical Society of the
	Pacific Conference Series, Vol. 302, Radio Pulsars, Bailes M., Nice
	D. J., Thorsett S. E., eds., p. 299
	
	[Kalogera et al(2004)Kalogera, Kim, Lorimer,
		Burgay, D'Amico, Possenti, Manchester, Lyne, Joshi, McLaughlin,
		Kramer, Sarkissian, & Camilo]Kalogera2004
	Kalogera V. et al., 2004, , 614, L137
	
	[Kasliwal et al(2017)Kasliwal, Nakar, Singer,
		Kaplan, Cook, Van Sistine, Lau, Fremling, Gottlieb, Jencson,
		Adams, Feindt, Hotokezaka, Ghosh, Perley, Yu, Piran, Allison,
		Anupama, Balasubramanian, Bannister, Bally, Barnes, Barway,
		Bellm, Bhalerao, Bhattacharya, Blagorodnova, Bloom, Brady,
		Cannella, Chatterjee, Cenko, Cobb, Copperwheat, Corsi, De,
		Dobie, Emery, Evans, Fox, Frail, Frohmaier, Goobar, Hallinan,
		Harrison, Helou, Hinderer, Ho, Horesh, Ip, Itoh, Kasen,
		Kim, Kuin, Kupfer, Lynch, Madsen, Mazzali, Miller, Mooley,
		Murphy, Ngeow, Nichols, Nissanke, Nugent, Ofek, Qi, Quimby,
		Rosswog, Rusu, Sadler, Schmidt, Sollerman, Steele, Williamson,
		Xu, Yan, Yatsu, Zhang, & Zhao]Kasliwal2017
	Kasliwal M. M. et al., 2017, Science, 358, 1559
	
	[Keane et al(2015)Keane, Bhattacharyya, Kramer,
		Stappers, Keane, Bhattacharyya, Kramer, Stappers, Bates,
		Burgay, Chatterjee, Champion, Eatough, Hessels, Janssen, Lee,
		van Leeuwen, Margueron, Oertel, Possenti, Ransom, Theureau, &
		Torne]Keane2015
	Keane E. et al., 2015, Advancing Astrophysics with the Square Kilometre
	Array (AASKA14), 40
	
	[Kim, Perera & McLaughlin(2015)Kim, Perera, &
		McLaughlin]Kim2015
	Kim C., Perera B. B. P., McLaughlin M. A., 2015, , 448, 928
	
	[Kiziltan & Thorsett(2010)]Kiziltan2010
	Kiziltan B., Thorsett S. E., 2010, , 715, 335
	
	[Kochanek & Piran(1993)]Kochanek1993
	Kochanek C. S., Piran T., 1993, , 417, L17
	
	[Komiya & Shigeyama(2016)]Komiya2016
	Komiya Y., Shigeyama T., 2016, , 830, 76
	
	[Kramer et al(2006)Kramer, Stairs, Manchester,
		McLaughlin, Lyne, Ferdman, Burgay, Lorimer, Possenti, D'Amico,
		Sarkissian, Hobbs, Reynolds, Freire, & Camilo]Kramer2006
	Kramer M. et al., 2006, Science, 314, 97
	
	[Lattimer & Schramm(1974)]Lattimer1974
	Lattimer J. M., Schramm D. N., 1974, , 192, L145
	
	[Lattimer & Schramm(1976)]Lattimer1976
	Lattimer J. M., Schramm D. N., 1976, , 210, 549
	
	[Lazarus et al(2016)Lazarus, Freire, Allen,
		Aulbert, Bock, Bogdanov, Brazier, Camilo, Cardoso, Chatterjee,
		Cordes, Crawford, Deneva, Eggenstein, Fehrmann, Ferdman,
		Hessels, Jenet, Karako-Argaman, Kaspi, Knispel, Lynch, van
			Leeuwen, Machenschalk, Madsen, McLaughlin, Patel, Ransom,
		Scholz, Seymour, Siemens, Spitler, Stairs, Stovall, Swiggum,
		Venkataraman, & Zhu]Lazarus2016
	Lazarus P. et al., 2016, , 831, 150
	
	[Leibler & Berger(2010)]Leibler2010
	Leibler C. N., Berger E., 2010, , 725, 1202
	
	[Lorimer et al(2006)Lorimer, Stairs, Freire,
		Cordes, Camilo, Faulkner, Lyne, Nice, Ransom, Arzoumanian,
		Manchester, Champion, van Leeuwen, Mclaughlin, Ramachandran,
		Hessels, Vlemmings, Deshpande, Bhat, Chatterjee, Han, Gaensler,
		Kasian, Deneva, Reid, Lazio, Kaspi, Crawford, Lommen, Backer,
		Kramer, Stappers, Hobbs, Possenti, D'Amico, &
		Burgay]Lorimer2006
	Lorimer D. R. et al., 2006, , 640, 428
	
	[Maoz, Mannucci & Nelemans(2014)Maoz, Mannucci, &
		Nelemans]Maoz2014
	Maoz D., Mannucci F., Nelemans G., 2014, , 52, 107
	
	[Martinez et al(2015)Martinez, Stovall, Freire,
		Deneva, Jenet, McLaughlin, Bagchi, Bates, &
		Ridolfi]Martinez2015
	Martinez J. G. et al., 2015, , 812, 143
	
	[Martinez et al(2017)Martinez, Stovall, Freire,
		Deneva, Tauris, Ridolfi, Wex, Jenet, McLaughlin, &
		Bagchi]Martinez2017
	Martinez J. G. et al., 2017, , 851, L29
	
	[Matsumoto, Nakar & Piran(2019)Matsumoto, Nakar, &
		Piran]Matsumoto2019
	Matsumoto T., Nakar E., Piran T., 2019, , 483, 1247
	
	[Matteucci et al(2014)Matteucci, Romano, Arcones,
		Korobkin, & Rosswog]Matteucci2014
	Matteucci F., Romano D., Arcones A., Korobkin O., Rosswog S., 2014,
	, 438, 2177
	
	[Osłowski et al(2011)Osłowski, Bulik,
		Gondek-Rosińska, & Belczyński]Oslowski2011
	Osłowski S., Bulik T., Gondek-Rosińska D., Belczyński K.,
	2011, , 413, 461
	
	[Piran(1992)]Piran1992
	Piran T., 1992, , 389, L45
	
	[Safarzadeh et al(2019)Safarzadeh, Ramirez-Ruiz,
		Andrews, Macias, Fragos, & Scannapieco]Safarzadeh2019
	Safarzadeh M., Ramirez-Ruiz E., Andrews J. J., Macias P., Fragos T.,
	Scannapieco E., 2019, , 872, 105
	
	[Shapiro & Teukolsky(1986)]ShapiroTeukolsky1986
	Shapiro S. L., Teukolsky S. A., 1986, Black Holes, White Dwarfs and
		Neutron Stars: The Physics of Compact Objects. p. 672
	
	[Simonetti et al(2019)Simonetti, Matteucci, Greggio,
		& Cescutti]Simonetti2019
	Simonetti P., Matteucci F., Greggio L., Cescutti G., 2019, arXiv
	e-prints
	
	[Snaith et al(2015)Snaith, Haywood, Di Matteo,
		Lehnert, Combes, Katz, & Gómez]Snaith2015
	Snaith O., Haywood M., Di Matteo P., Lehnert M. D., Combes F., Katz
	D., Gómez A., 2015, , 578, A87
	
	[Stairs et al(2002)Stairs, Thorsett, Taylor, &
		Wolszczan]Stairs2002
	Stairs I. H., Thorsett S. E., Taylor J. H., Wolszczan A., 2002, ,
	581, 501
	
	[Stovall et al(2018)Stovall, Freire, Chatterjee,
		Demorest, Lorimer, McLaughlin, Pol, van Leeuwen, Wharton,
		Allen, Boyce, Brazier, Caballero, Camilo, Camuccio, Cordes,
		Crawford, Deneva, Ferdman, Hessels, Jenet, Kaspi, Knispel,
		Lazarus, Lynch, Parent, Patel, Pleunis, Ransom, Scholz,
		Seymour, Siemens, Stairs, Swiggum, & Zhu]Stovall2018
	Stovall K. et al., 2018, , 854, L22
	
	[Suda et al(2008)Suda, Katsuta, Yamada, Suwa,
		Ishizuka, Komiya, Sorai, Aikawa, & Fujimoto]Suda2008
	Suda T. et al., 2008, , 60, 1159
	
	[Swiggum et al(2015)Swiggum, Rosen, McLaughlin,
		Lorimer, Heatherly, Lynch, Scoles, Hockett, Filik, Marlowe,
		Barlow, Weaver, Hilzendeger, Ernst, Crowley, Stone, Miller,
		Nunez, Trevino, Doehler, Cramer, Yencsik, Thorley, Andrews,
		Laws, Wenger, Teter, Snyder, Dittmann, Gray, Carter, McGough,
		Dydiw, Pruett, Fink, & Vanderhout]Swiggum2015
	Swiggum J. K. et al., 2015, , 805, 156
	
	[Tauris et al(2017)Tauris, Kramer, Freire, Wex,
		Janka, Langer, Podsiadlowski, Bozzo, Chaty, Kruckow, van den
			Heuvel, Antoniadis, Breton, & Champion]Tauris2017
	Tauris T. M. et al., 2017, , 846, 170
	
	[Totani et al(2008)Totani, Morokuma, Oda, Doi, &
		Yasuda]Totani2008
	Totani T., Morokuma T., Oda T., Doi M., Yasuda N., 2008, , 60,
	1327
	
	[Tsujimoto & Shigeyama(2014)]Tsujimoto2014
	Tsujimoto T., Shigeyama T., 2014, , 565, L5
	
	[van Leeuwen et al(2015)van Leeuwen, Kasian, Stairs,
		Lorimer, Camilo, Chatterjee, Cognard, Desvignes, Freire,
		Janssen, Kramer, Lyne, Nice, Ransom, Stappers, &
		Weisberg]vanLeeuwen2015
	van Leeuwen J. et al., 2015, , 798, 118
	
	[Vuong(1989)]Vuong89
	Vuong Q. H., 1989, Econometrica, 57, 307
	
	[Wallner et al(2015)Wallner, Faestermann, Feige,
		Feldstein, Knie, Korschinek, Kutschera, Ofan, Paul, Quinto,
		Rugel, & Steier]Wallner2015
	Wallner A. et al., 2015, Nature Communications, 6, 5956
	
	[Wanderman & Piran(2015)]WP2015
	Wanderman D., Piran T., 2015, , 448, 3026
	
	[Wehmeyer, Pignatari & Thielemann(2015)Wehmeyer, Pignatari,
		& Thielemann]Wehmeyer2015
	Wehmeyer B., Pignatari M., Thielemann F.-K., 2015, , 452, 1970
	
	[Weisberg, Nice & Taylor(2010)Weisberg, Nice, &
		Taylor]Weisberg2010
	Weisberg J. M., Nice D. J., Taylor J. H., 2010, , 722, 1030
	
	[Weisz et al(2015)Weisz, Dolphin, Skillman,
		Holtzman, Gilbert, Dalcanton, & Williams]Weisz2015
	Weisz D. R., Dolphin A. E., Skillman E. D., Holtzman J., Gilbert
	K. M., Dalcanton J. J., Williams B. F., 2015, , 804, 136
	
	[Yungelson & Livio(2000)]Yungelson2000
	Yungelson L. R., Livio M., 2000, , 528, 108
	


	
	
	
	
	
	
