Measuring the delay time distribution of binary neutron stars. I. Through Scaling Relations of the Host Galaxies of Gravitational Wave Events

TL;DR

This paper proposes a method to constrain the delay time distribution of binary neutron star mergers using host galaxy demographics from gravitational wave detections, requiring about a thousand galaxy observations.

Contribution

It introduces a new approach leveraging galaxy mass and star formation history correlations to estimate DTD parameters from GW host galaxy data.

Findings

Approximately 1,000 host galaxies are needed to constrain the DTD.

Massive galaxies primarily encode the DTD information.

Fewer galaxies (~100) suffice if the minimum delay time is fixed.

Abstract

The delay time distribution of (DTD) of binary neutron stars (BNS) remains poorly constrained, mainly by the small known population of Galactic binaries, the properties of short gamma-ray burst host galaxies, and inferences from $r$-process enrichment. In the new era of BNS merger detections through gravitational waves (GW), a new route to the DTD is the demographics of the host galaxies, localized through associated electromagnetic counterparts. This approach takes advantage of the correlation between star formation history (SFH) and galaxy mass, such that the convolution of the SFH and DTD impacts the BNS merger rate as a function of galaxy mass. Here we quantify this approach for a power law DTD governed by two parameters: the power law index ($\Gamma$) and a minimum delay time ($t_{\rm min}$). Under the reasonable assumption that EM counterparts are likely only detectable in the local universe, accessible by the current generation of GW detectors, we study how many host galaxies at $z\sim 0$ are required to constrain the DTD parameters. We find that the DTD is mainly imprinted in the statistics of massive galaxies (stellar mass of $M_*\gtrsim 10^{10.5}$ M$_\odot$, comparable to the host galaxy of GW170817). Taking account of relevant uncertainties we find that $\mathcal{O}(10^3)$ host galaxies are required to constrain the DTD; for a fixed value of $t_{\rm min}$, as done in previous analyses of the DTD, $\mathcal{O}(10^2)$ host galaxies will suffice. Such a sample might become available within the next two decades, prior to the advent of third-generation GW detectors.

Equations

$$\dot{n}(M_h)=$$

$$\frac{dt}{dz}=\frac{-1}{(1+z)E(z)H_0},$$

$$E(z)=\sqrt{{\Omega }_{m,0}(1+z{)}^3+{\Omega }_{k,0}(1+z{)}^2+{\Omega }_{\Lambda }(z)}.$$

$$\varphi (M_h)\equiv \frac{dn}{d{M}_h}=\frac{\bar{\rho }}{M_h}f(\nu )\frac{d\nu }{d{M}_h},$$

$$\nu f(\nu )=A\sqrt{\frac{a\nu }{2\pi }}\left[1+(a\nu {)}^{-p}\right]e^{-a\nu /2},$$

$${\sigma }_{\chi }^2(M,0)={\sigma }_{\chi }^2(R,0)={\int }_0^{\infty }\frac{dk}{2{\pi }^2}{k}^2{P}_{\chi }(k)w^2(kR)$$

$$D(z)\equiv\frac{H(z)}{H(0)}\int_{z}^{\infty}\frac{dz^{\prime}(1+z^{\prime})}{H^{3}(z^{\prime})}\Bigg{[}\int_{0}^{\infty}\frac{dz^{\prime}(1+z^{\prime})}{H^{3}(z^{\prime})}\Bigg{]}^{-1}.$$

$$f(>M_h)=\frac{{\int }_{M_h}^{M_{h,\max }}\varphi (M_h^{\prime })\dot{n}(M_h^{\prime })d{M}_h^{\prime }}{{\int }_{M_{h,\min }}^{M_{h,\max }}\varphi (M_h^{\prime })\dot{n}(M_h^{\prime })d{M}_h^{\prime }},$$

Full text

1 Introduction

The delay time distribution of binary neutron stars is currently only weakly constrained, mainly by the statistics of the small known sample of Galactic binary neutron stars (e.g., Vigna-Gómez et al. 2018), from arguments related to $r$-process enrichment (e.g., Matteucci et al., 2014; Safarzadeh et al., 2018), and from the properties of short gamma-ray burst (SGRB) host galaxies (Zheng & Ramirez-Ruiz, 2007; Leibler & Berger, 2010; Fong et al., 2013; Berger, 2014). The DTD is generally expected to follow a power law distribution based on the following reasoning. After the formation of the BNS, the binary’s orbit decays through the emission of gravitational waves on a timescale that depends on the binary’s separation as $t\propto a^{4}$, where $a$ is the semi-major axis of the binary at formation (Peters, 1964). The distribution of the merging times, therefore, depends on the distribution of the initial orbital separation modeled as $dN/da\propto a^{-\beta}$. The initial distribution of the O/B star (the progenitors of the NSs) is assumed to follow a power law $dN/da\propto a^{-1}$. If the binary goes through a common envelope phase, then the distribution of the separation becomes steeper and approaches $dN/da\propto a^{-3}$. Therefore, the expected merger times follow $dN/dt_{\rm merge}\propto t^{\Gamma}$, where $\Gamma\equiv-\beta/4-3/4$. For those two limiting cases, $\Gamma$ ranges from $-1.5$ to $-1$ (e.g., Belczynski et al., 2018). The various weak observational constraints are roughly in agreement with these values.

Separate from the slope of the power law distribution, the minimum timescale for BNS mergers ($t_{\rm min}$) is another parameter that is equally important in determining the merger rate across cosmic time. From population synthesis models $t_{\rm min}$ could be as short as a few Myr (Dominik et al., 2012), but various effects could serve to set a minimum initial separation that will increase the value of $t_{\rm min}$. Observationally, the two key DTD parameters are approximately degenerate with each other in that it is not trivial to distinguish between a DTD with a steep slope but larger $t_{\rm min}$ and a DTD with shallow slope but shorter $t_{\rm min}$. Recent simulations have shown that fast merging channels are needed to explain the fraction of all the metal-poor stars that are $r$-process enriched (Matteucci et al., 2014; Safarzadeh et al., 2018), as well as for $r$-process enrichment of ultra-faint dwarf galaxies (Safarzadeh & Scannapieco, 2017; Safarzadeh et al., 2019). We note that most previous attempts to determine the DTD have assumed a small value of $t_{\rm min}$, leaving $\Gamma$ as the only free parameter. Such an assumption is also made in the case of the DTD of Type Ia supernovae (e.g., $t_{\rm min}\sim 40$ Myr based on the minimum lifetime of stars that produce white dwarfs; Maoz et al. 2012, 2014), leading to better constraints on the power law index.

The shape of the DTD is also imprinted in the demographics of the galaxies that host BNS mergers in the local universe, manifested either as the host galaxies of SGRBs (Berger, 2014), or as the host galaxies of gravitational wave events that can in turn be pinpointed through the detection of associated electromagnetic counterparts (e.g., kilonovae; Abbott et al. 2017; Coulter et al. 2017; Soares-Santos et al. 2017). This is because the star formation histories of galaxies are determined by their masses, and the convolution of the SFH with the DTD will therefore impact the mass distribution of BNS merger host galaxies (Zheng & Ramirez-Ruiz, 2007; Leibler & Berger, 2010; Fong et al., 2013; Artale et al., 2019). The detection of the BNS merger GW170817, and the identification of its host galaxy, pave the way for utilizing this approach to constrain the DTD.

Here we use galaxy scaling relations to explore the impact of the DTD on the distribution of BNS merger host galaxies, and explore the number of events required to constrain the DTD. The upcoming observing campaigns with Advanced LIGO/Virgo and the upcoming detectors KAGRA and IndIGO are expected to yield BNS merger samples of $\mathcal{O}(10)$, $\mathcal{O}(100)$, and $\mathcal{O}(1000)$ within the next year, $\sim 5$ years, and $\sim 20$ years, respectively, before the advent of third-generation GW detectors. The structure of this Letter is as follows: In §2 we demonstrate how the shape of the DTD affects the demographics of BNS merger host galaxies in the local universe; in §3 we explore the sample size required to constrain the shape of the DTD; and in §4 we discuss the caveats involved in this analysis. We adopt the Planck 2015 cosmological parameters (Planck Collaboration et al., 2016): $\Omega_{M}=0.308$, $\Omega_{\Lambda}=0.692$, $\Omega_{b}=0.048$, and $H_{0}=0.678$ km s*-1* Mpc*-1*.

2 Method

We can write the BNS merger rate for a galaxy with halo mass, $M_{h}$, at $z=0$ as:

[TABLE]

where

[TABLE]

and

[TABLE]

Here, $\psi(M_{h},z)$ is the mean star formation history (SFH) of a galaxy with $M_{h}$ at $z=0$ parametrized following Moster et al. (2012). We integrate the SFH from $z_{b}=10$ to [math] (where the choice of maximum redshift has little impact on the calculation); $t_{b}$ is the cosmic time corresponding to $z_{b}$; $\lambda$ is the BNS mass efficiency, assumed to be a fixed value of $10^{-5}M_{\odot}^{-1}$ independent of redshift or environment; $dP_{m}/dt$ is the merger time distribution, which we parametrize to follow a power law, $\propto t^{\Gamma}$ with a minimum delay time, $t_{\rm min}$ and not evolving with redshift. Although the DTD for binary black holes is likely highly dependent on metallicity, the DTD for BNS systems has been argued to be at most weakly dependent on metallicity (Dominik et al., 2012). The mass efficiency is assumed to be constant, although this parameter could be fit for in principle (Safarzadeh et al., 2019), due to the limited depth of adLIGO it acts as a normalization constant that could be ignored when studying the distribution of host galaxy masses in the local universe. We note that the delay time refers to the time since birth of the ZAMS stars and not when the BNS is formed. We also impose a maximum delay time of 10 Gyr for our fiducial case, but our results are not sensitive for a longer maximum delay time.

We compute the merger rate per galaxy as a function of halo mass for a grid of 9 joint choices of $\Gamma=[-1.5,-1.0,-0.5]$ and $t_{\rm min}=[10,100,1000]$ Myr. In Figure 1 we show the predicted merger rate ($\dot{n}$) as a function of halo mass and stellar mass for the 9 DTDs. We find that the key difference between the various DTDs is apparent at $M_{h}\gtrsim 10^{12}$ M⊙, corresponding to $M_{*}\gtrsim 10^{10.5}$ M⊙. This is primarily because on average galaxies of lower masses have fairly flat star formation histories that are therefore not sensitive to convolution with the different DTDs. The high mass galaxies, on the other hand, have star formation histories that peak at progressively earlier cosmic time with larger mass. Therefore, we find that DTDs that favor long merger timescales (e.g., $\Gamma=-1/2$ and $t_{\rm min}=1$ Gyr) lead to a higher representation of massive host galaxies.

To determine the observed mass distribution of BNS merger host galaxies we need to rescale $\dot{n}$ with the halo mass function (HMF), computed following Press & Schechter (1974):

[TABLE]

where $n$ is the number density of haloes, $\nu$ is the peak-height of perturbations, $\bar{\rho}$ is the average density of the universe, and the first crossing distribution, $f(\nu)$ (Bond et al., 1991), is obtained from the ellipsoidal collapse model as:

[TABLE]

with $A=0.322$, $p=0.3$, and $a=0.75$ (Sheth & Tormen, 2002). Here, the peak height, $\nu$, is defined as $\nu\equiv\delta_{c,0}^{2}{\sigma_{\chi}(R,z)}^{-2},$ with $\delta_{c,0}=1.686$. The variance is $\sigma^{2}_{\chi}(M,z)=\sigma^{2}_{\chi}(M,0)D(z)^{2},$ with

[TABLE]

where $M=4\pi R^{3}\Omega_{M}\rho_{c}/3$, $w(kR)\equiv 3j_{1}(kR)/(kR),$ with $j_{1}(x)\equiv(\sin x-x\cos x)/x^{2},$ and $D(z)$ is the linear growth factor

[TABLE]

We compute the cumulative fraction of the BNS merger host halos with mass above $M_{h}$ as:

[TABLE]

where we consider halos in the range $M_{h,\rm min}=10^{11}$ to $M_{h,\rm max}=10^{14}$ M⊙. In Figure 1 we show the cumulative distribution function (CDF) for BNS merger host galaxies as a function of $M_{h}$ and $M_{*}$ for the 9 DTDs. We find that there is a difference of about 0.5 dex in the median value of $M_{h}$ for the range of DTDs, and about 0.7 dex in the value of $M_{h}$ for the top 20% most massive galaxies; at the low mass end the CDFs converge. Similarly, in terms of stellar mass we find a nearly order of magnitude spread in the median value of $M_{*}$. In both CDFs we again find a clear degeneracy between $\Gamma$ and $t_{\rm min}$, such that a DTD with a shallower power law index and small value of $t_{\rm min}$ is similar in shape to one with a steeper power law index and a large $t_{\rm min}$.

Finally, to determine the sample size needed to constrain the DTD we draw from the constructed CDF of a given pair of [$\Gamma$, $t_{\rm min}$], and perform a Kolmogorov-Smirnov (KS) test against other possible CDFs constructed in a $10\times 10$ interpolated 2D plane of $\Gamma-t_{\rm min}$. We use a threshold $P<0.01$ to consider the CDFs as being drawn from different underlying distributions. For each case we repeat the KS test 10 times and determine the median value of $P$.

We note that the overall approach requires an identified EM counterpart and host galaxy; with only two examples to date (GW170817 and likely S190425z; Hosseinzadeh et al. 2019) it is difficult to assess the counterpart identification success rate. We show in Safarzadeh et al. (2019) how future third-generation GW detectors could enable a measurement of the DTD without the need for EM counterparts. Moreover, we focused on BNS systems formed in the field, while there are other proposed mechanism for their formation, although those are expected to be a minor contributor (e.g., Grindlay et al. 2006; Lee et al. 2010).

3 Results

In Figure 2 we plot the value of $\log P$ for three different injected DTDs and three different sample sizes (30, 100, and 1000 BNS merger host galaxies) in the plane of $\Gamma-t_{\rm min}$. This figure uses the halo mass CDFs shown in Figure 1. We find that a sample size of $\mathcal{O}(10^{2})$ is required to begin to rule out portions of the DTD parameter space. With a sample size of $\mathcal{O}(10^{3})$ a significant portion of the parameter space can be ruled out, but with a possible remaining degeneracy between $\Gamma$ and $t_{\rm min}$ depending on the true values of the parameters.

So far we have cast our results in terms of $M_{h}$, but observationally we determine the stellar mass ($M_{*}$) based on modeling a galaxy’s spectral energy distribution. While there is uncertainty involved in estimating $M_{*}$, the connection to $M_{h}$, based on abundance matching techniques (Behroozi et al., 2013) is likely to dominate. For example, Blanchard et al. (2017) determined the stellar mass of the host galaxy of GW170817 (NGC 4993) to be $\log(M_{*}/M_{\odot})=10.90^{+0.03}_{-0.03}$ by reconstructing its star formation history using UV to IR data. However, the corresponding halo mass is estimated to be $\log(M_{h}/M_{\odot})=13.0^{+0.3}_{-0.3}$, with the much larger uncertainty due to the dispersion in the $M_{h}-M_{*}$ relation (Behroozi et al., 2013). To accommodate this uncertainty into our analysis, we disperse the halo masses by 0.3 dex when sampling from a given CDF, and then follow the same procedure as above to determine $P$ values. The results of including this systematic uncertainty of 0.3 dex are equivalent to the uncertainty on the SFH of an observed galaxy with a given stellar mass. The inclusion of this additional uncertainty reduces the constraining power of the observed host galaxy sample on the DTD. However, in Figure 3 we show that a sample size of $\mathcal{O}(10^{3})$ can still provide significant constraints on the DTD even when this systematic uncertainty is included.

The overall statistics of the fraction of the $\Gamma-t_{\rm min}$ parameter space that can be ruled out with different sample sizes is shown in Figure 4. Here we inject 50 random [$\Gamma,t_{\rm min}$] pairs, perform the same KS test analysis and determine the fraction of parameter space with $P<0.01$ (i.e., the excluded fraction). We find that the median excluded fraction of the parameter space is about 10%, 25%, and 65% for sample sizes of 30, 100, and 1000 BNS merger host galaxies, respectively.

We note that our analysis uses a DTD with two free parameters, while generally only $\Gamma$ is used as a free parameter, and $t_{\rm min}$ is fixed at a small value ($\sim 10$ Myr). This simplifying assumption is based on the notion that some binaries can merge as soon as the second neutron star is formed. This approach has also been used in analyses of the DTD of Type Ia supernovae, in which a value of $t_{\rm min}\sim 40$ Myr is often assumed (e.g., Maoz et al. 2012, 2014), leaving only $\Gamma$ as a free parameter. In Figure 5 we show the resulting constraints on $\Gamma$ if we fix $t_{\rm min}=10$ Myr and repeat our analysis. We find that in this simplified model the value of $\Gamma$ can be determined with $\approx 30\%$ uncertainty with a reduced sample size of $\sim 300$ BNS merger host galaxies. Such a sample can be accumulated in about one-third of the time compared to the requirement when both $\Gamma$ and $t_{\rm min}$ are free parameters.

4 Summary

We showed how the DTD of BNS systems can be constrained through the demographics of the host galaxies of BNS mergers detected in gravitational waves and pinpointed through electromagnetic observations. We focused on the case of a DTD parameterized as a power law, although other possible DTD shapes have been proposed in the literature (e.g., Simonetti et al. 2019). Our analysis is similar to that of Zheng & Ramirez-Ruiz (2007), proposed in the context of SGRB host galaxy demographics, but with the difference that we model the DTD with two parameters, while those authors used just $\Gamma$ and fixed $t_{\rm min}$.

Our results show that $\mathcal{O}(10^{3})$ host galaxies are needed to constrain a two-parameter DTD, although some degeneracy between $\Gamma$ and $t_{\rm min}$ is intrinsically difficult to resolve even with a large sample of events. On average, about two-thirds of the $\Gamma-t_{\rm min}$ parameter space can be ruled out with such a sample. In the case when only $\Gamma$ is a free parameter, a sample size of about 300 BNS merger host galaxies is sufficient for a 30% uncertainty on $\Gamma$. The current range of BNS merger rates from Advanced LIGO/Virgo Observing Runs 1 and 2 is $110-3840$ Gpc*-3* yr*-1* (The LIGO Scientific Collaboration et al., 2018), indicating that a sample of $\mathcal{O}(10^{3})$ events might be achieved within a couple of decades at a design sensitivity of 200 Mpc, or potentially even faster in the case of A+ (Barsotti et al., 2018), with an expected factor of 2 increase in the BNS merger detection range. Thus, a constraint on the DTD using the demographics of BNS merger host galaxies can be achieved before the advent of the third-generation GW detectors. We note that a similar approach using SGRB host galaxies is likely to take longer given an identification rate of only a few events per year (Berger, 2014).

In an upcoming paper we will investigate a related approach to constraining the DTD, using the measured individual SFH of BNS merger host galaxies (rather than the mean scaling relations used here). This is similar to the approach used by Blanchard et al. (2017) for the host galaxy of GW170817, and by Maoz & Graur (2017) in the context of Type Ia SNe.

We are thankful to Enrico Ramirez-Ruiz, Or Graur, and Evan Scannapieco for helpful discussions. This work was supported by the National Science Foundation under grant AST14-07835 and by NASA under theory grant NNX15AK82G. The Berger Time-Domain Group at Harvard is supported in part by NSF under grant AST-1714498 and by NASA under grant NNX15AE50G. MTS is thankful to Harvard-Smithsonian Center for Astrophysics for hospitality which made this work possible.