Measuring the Delay Time Distribution of Binary Neutron Stars. II. Using the Redshift Distribution from Third-generation Gravitational-wave Detectors Network

The joint gravitational wave (GW) and electromagnetic (EM) detections of the binary neutron star (BNS) merger, GW170817 (Abbott et al. 2017b), marked the dawn of multi-messenger astronomy. As the current generation of GW detectors increases in sensitivity and number, the local rate of binary neutron star (BNS) mergers will soon be determined accurately for the first time, providing initial insight into the formation channels of these binaries. Currently, the BNS merger rate is weakly constrained by the single detection of GW170817 (${1540}_{-1220}^{+3200}$ Gpc⁻³ yr⁻¹; Abbott et al. 2017c), by the small known sample of Galactic BNS systems (${21}_{-14}^{+28}$ Myr⁻¹; Kim et al. 2015), and by the beaming-corrected rate of short gamma-ray bursts (SGRBs; ${270}_{-180}^{+1580}$ Gpc⁻³ yr⁻¹; Fong et al. 2015). These rates are in broad agreement, but the uncertainties from all methods span at least two orders of magnitude.

Still, even when the local BNS merger rate is well determined, the more fundamental distribution of merger delay times may not be. The delay time distribution (DTD) encodes the time span between the formation of the BNS system (or alternatively the time since the formation of the parent stars) until the two neutron stars merge through the emission of GWs. The DTD therefore provides fundamental insight into the evolutionary processes that govern the initial separation of the binaries, including poorly understood effects such as common envelope evolution.

The DTD is usually parametrized as a power-law distribution above some minimum merger timescale, t_min, based on the following arguments: after the BNS formation, the binary's orbit decays due to the emission of GWs on a timescale that depends on the initial semimajor axis (a) as t ∝ a⁴. Therefore, the resulting distribution of the merger times depends on the distribution of initial semimajor axes, dN/da ∝ a^−β. The initial semimajor axis distribution of the O/B stellar progenitors is assumed to follow a power law dN/da ∝ a⁻¹. If the binary experiences a common envelope phase, then the distribution becomes steeper. Therefore, the expected merger times follow dN/dt_merge ∝ t^{−β/4−3/4}, where we define Γ ≡ −β/4 − 3/4 (Belczynski et al. 2018).

Insight on the form of the DTD has been gained from studies of the small population of Galactic BNS systems (Vigna-Gómez et al. 2018), from the properties of SGRB host galaxies (Zheng & Ramirez-Ruiz 2007; O'Shaughnessy et al. 2008; Leibler & Berger 2010; Fong et al. 2013; Behroozi et al. 2014; Berger 2014), and from arguments related to r-process enrichment (Komiya et al. 2014; Matteucci et al. 2014; Shen et al. 2015; van de Voort et al. 2015; Côté et al. 2018; Hotokezaka et al. 2018; Safarzadeh et al. 2018, 2019). These results point to the need for a fast-merging channel if BNSs are assumed to be the primary source of r-process enrichment in the universe, which suggest that t_min may be rather small, ≲0.1 Gyr, or the slope of the DTD is steep. Population synthesis models have also made various predictions for the values of Γ and t_min, but those are dependent on uncertain binary evolution processes (Dominik et al. 2012). We stress that the local BNS merger rate in itself cannot fully characterize the DTD because it also depends on an additional unknown parameter: the efficiency of BNS formation per unit stellar mass, λ.

In a recent Letter, Safarzadeh & Berger (2019; hereafter, Paper I) showed that the mass distribution of BNS merger host galaxies at z ≈ 0 can provide insight into Γ and t_min with a sample size of ${ \mathcal O }({10}^{2}-{10}^{3})$. This is based on the fact that, on average, galaxy star formation histories (SFHs) depend on their mass, and hence the convolution of the DTD and SFH leads to a specific prediction about the mass function of BNS merger host galaxies. Such an observational approach to determining the DTD is only feasible in the local universe due to the required detection of EM counterparts that will in turn lead to the identification of the host galaxies. It is anticipated that Advanced Laser Interferometer Gravitational-Wave Observatory (LIGO)/Virgo, joined by the Kamioka Gravitational Wave Detector (KAGRA)⁶ and the Indian Initiative in Gravitational-wave Observations (IndIGO),⁷ can produce the required sample size within the next two decades.

Here, we instead explore how the DTD can be determined by directly observing the redshift distribution of BNS mergers well beyond the local universe. Mapping the rate of BNS mergers as a function of redshift can break the degeneracy between the shape of the DTD (Γ and t_min) and the BNS mass efficiency (λ) when comparing to the cosmic SFH. This approach requires an order of magnitude increase in GW detector sensitivity to detect BNS mergers at cosmological distances. Such an improvement is expected for third-generation ground-based observatories such as Einstein Telescope⁸ (ET; Punturo et al. 2010) and Cosmic Explorer⁹ (CE; Abbott et al. 2017a). However, at these distances, it is unlikely that EM counterparts will be detected for the majority of events, and therefore the distance (redshift) information will rely directly on the GW signal itself. We explore how the inherent distance-inclination degeneracy affects the ability to determine the DTD.

The structure of this Letter is as follows. In Section 2 we delineate the method of estimating the observed redshift distribution of BNS mergers as a function of DTD and mass efficiency, for different GW interferometer networks; in Section 3 we show the results of DTD determination for existing GW detectors, which are only sensitive to the local universe; in Section 4 we expand our analysis to a future network of ET and CE, including a determination of the expected redshift uncertainties, and show the resulting constraints on the DTD and mass efficiency. We discuss some caveats and summarize the key results in Section 5. We adopt the Planck 2015 cosmological parameters (Planck Collaboration et al. 2016) where Ω_M = 0.308, Ω_Λ = 0.692, Ω_b = 0.048 are total matter, vacuum, and baryonic densities, in units of the critical density, ρ_c, H₀ = 67.8 km s⁻¹ Mpc⁻¹ is the Hubble constant, and σ₈ = 0.82 is the variance of linear fluctuations on the 8 h⁻¹ Mpc scale.

The BNS merger rate as a function of redshift is a convolution of the DTD with the cosmic star formation rate density (SFRD):

$$\begin{eqnarray}&&\dot{n}(z)={\int }_{{z}_{b}=10}^{{z}_{b}=z}\lambda \displaystyle \frac{{{dP}}_{m}}{{dt}}(t-{t}_{b}-{t}_{\min })\psi ({z}_{b})\displaystyle \frac{{dt}}{{dz}}({z}_{b}){{dz}}_{b},\end{eqnarray} \tag{ 1 }$$

where dt/dz = −[(1 + z)E(z)H₀]⁻¹, and $E(z)\,=\sqrt{{{\rm{\Omega }}}_{m,0}{\left(1+z\right)}^{3}+{{\rm{\Omega }}}_{k,0}{\left(1+z\right)}^{2}+{{\rm{\Omega }}}_{{\rm{\Lambda }}}(z)}$. Here, λ is the currently unknown BNS mass efficiency (assumed not to evolve¹⁰ with redshift) used as a free parameter that we try to recover alongside the parameters governing the DTD; t_b is the time corresponding to the redshift z_b; dP_m/dt is the DTD, parametrized to follow a power-law distribution (∝t^Γ) with a minimum delay time, t_min that refers to the time since birth of the zero-age main-sequence stars and not when the BNS system formed. Therefore, t_min corresponds to the sum of the nuclear lifetime of the lowest-mass component of the binary system and the minimal gravitational delay that is induced by the existence of a minimal separation between the two newly born neutron stars. We also impose a maximum delay time of 10 Gyr for our fiducial case, although this does not affect our results, although we note that more than half of the observed BNS systems in the Milky Way (MW) have merger times that are more than 10 Gyr (Pol et al. 2019). We adopt the cosmic SFRD¹¹ from Madau & Dickinson (2014):

$$\begin{eqnarray}&&\psi (z)=0.015\displaystyle \frac{{\left(1+z\right)}^{2.7}}{1+{\left[(1+z)/2.9\right]}^{5.6}}\,{M}_{\odot }{\mathrm{yr}}^{-1}\,{\mathrm{Mpc}}^{-3}.\end{eqnarray} \tag{ 2 }$$

To determine the observed BNS merger rate as a function of redshift we need to consider the matched filtering signal-to-noise ratio (S/N) as a function of GW detector sensitivity (Finn 1996):

$$\begin{eqnarray}&&\rho (z)=8{\rm{\Theta }}\displaystyle \frac{{r}_{0}}{{D}_{L}}{\left(\displaystyle \frac{{{ \mathcal M }}_{z}}{1.2{M}_{\odot }}\right)}^{5/6}\sqrt{\zeta ({f}_{\max })},\end{eqnarray} \tag{ 3 }$$

where ${{ \mathcal M }}_{z}=(1+z){ \mathcal M }$ is the redshifted chirp mass, D_L is the luminosity distance, Θ is the orientation function, and

$$\begin{eqnarray}\begin{array}{rcl}{r}_{0}^{2} & \equiv & \displaystyle \frac{5}{192\pi }{\left(\displaystyle \frac{3G}{20}\right)}^{5/3}{x}_{7/3}\displaystyle \frac{{M}_{\odot }^{2}}{{c}^{3}},\\ {x}_{7/3} & \equiv & {\displaystyle \int }_{0}^{\infty }\displaystyle \frac{{df}{\left(\pi {M}_{\odot }\right)}^{2}}{{\left(\pi {{fM}}_{\odot }\right)}^{7/3}{S}_{h}(f)},\\ \zeta ({f}_{\max }) & \equiv & \displaystyle \frac{1}{{x}_{7/3}}{\displaystyle \int }_{0}^{2{f}_{\max }}\displaystyle \frac{{df}{\left(\pi {M}_{\odot }\right)}^{2}}{{\left(\pi {{fM}}_{\odot }\right)}^{7/3}{S}_{h}(f)},\end{array}\end{eqnarray} \tag{ 4 }$$

where 2f_max is the wave frequency at which the inspiral detection template ends, r₀ denotes the characteristic distance sensitivity, and S_h(f) is the detector's noise power spectral density. The intrinsic chirp mass, ${ \mathcal M }$, is given in terms of the component masses by

$$\begin{eqnarray}&&{ \mathcal M }={\left(\displaystyle \frac{{m}_{1}{m}_{2}}{{({m}_{1}+{m}_{2})}^{2}}\right)}^{3/5}({m}_{1}+{m}_{2}).\end{eqnarray} \tag{ 5 }$$

Here we assume that both neutron stars have mass of m₁ = m₂ = 1.4 M_⊙. The frequency at the end of the inspiral (taken to correspond to the innermost stable circular orbit) is

$$\begin{eqnarray}&&{f}_{\max }=\displaystyle \frac{785\ \mathrm{Hz}}{1\,+\,z}\left(\displaystyle \frac{2.8\,{M}_{\odot }}{M}\right),\end{eqnarray} \tag{ 6 }$$

where M is the total mass of the binary. In Figure 1 we show the sensitivity curves for Advanced LIGO (AdLIGO), ET, and CE. The substantial reduction in noise amplitude for the third-generation detectors with respect to AdLIGO leads to an increase in the typical values of r₀ from ≈0.1 to ≈1.5 Gpc. Finally, the observed BNS merger rate as a function of redshift is given by

$$\begin{eqnarray}&&{R}_{D}(z)=\displaystyle \frac{{{dV}}_{c}}{{dz}}\displaystyle \frac{\dot{n}(z)}{1\,+\,z}{P}_{\det }(z),\end{eqnarray} \tag{ 7 }$$

where P_det(z) is defined in Appendix A, and the redshift derivative of the comoving volume is given by ${{dV}}_{c}/{dz}\,=(4\pi c/{H}_{0})[{D}_{L}^{2}/{\left(1+z\right)}^{2}E(z)]$.

Zoom In Zoom Out Reset image size

In Figure 2 we show the intrinsic merger rate density, $\dot{n}$(z), for nine different choices of the DTD, with Γ = [−1.5, −1, −0.5] and t_min = [10,100,1000] Myr, and a fixed mass-efficiency value of λ = 10⁻⁵ ${M}_{\odot }^{-1}$. For comparison, we also show the curve corresponding to no delay (i.e., the cosmic SFRD). Clearly, DTDs that prefer longer delays result in a merger distribution that is skewed to lower redshifts, with a higher merger rate at z ≈ 0, but with some degeneracy between Γ and t_min. However, because the value of λ is not presently known, all of the DTDs can reproduce the same local rate by simply scaling λ appropriately. This is essentially why a local measurement of the merger rate cannot by itself constrain the DTD.

Zoom In Zoom Out Reset image size

To explore how well current and third-generation GW detectors can determine the DTD, we inject a specific DTD model (Γ, t_min, λ), generate the resulting redshift distribution with associated uncertainties, and then fit this distribution using an interpolation table that is based on the nine input DTDs. We fit for the input parameters using Markov Chain Monte Carlo (MCMC) sampling with emcee, a python-based affine invariant sampler (Foreman-Mackey et al. 2013). The likelihood function is ln(L) = −χ²/2, with ${\chi }^{2}={\sum }_{i=0}^{i=N}({R}_{D,i}-\hat{{R}_{D,i}})/{\sigma }_{t,i}^{2}$. Here the summation is over all of the redshift bins; R_D,i and $\hat{{R}_{D,i}}$ are the constructed and simulated detection rates at redshift bin i, respectively; σ_t,i is the total error on the detection rate at redshift bin i, which is a combination of the Poisson error and the error due to the distance-inclination degeneracy from GW data, ${\sigma }_{t}^{2}={\sigma }_{p}^{2}+{\sigma }_{z}^{2};$ ${\sigma }_{p}=\sqrt{N}$, where N is the expected number of events at a given redshift during the integrated observation time of length T_obs; and the redshift uncertainty (σ_z) for each redshift bin is estimated based on the vertical distance from the mean expected detection rate to the upper envelope corresponding to when the detections' redshift are all biased high. We model the redshift uncertainty as δz/z = 0.1z based on rescaled simulations of BBH redshift uncertainty estimates as detailed in Appendix B. We adopt a flat prior distribution for all of our parameters in the logλ ∈ [−7, −3], log t_min ∈ [1,3], and Γ ∈ [−1.5, −0.5].

The predicted observed redshift distribution for the current generation of detectors at design sensitivity is shown in Figure 3. As expected, because the detection distance is limited to only a few hundred Mpc, all of the DTDs predict the same shape of observed distribution, with a simple change in scaling that can be accommodated by varying the unknown value of λ.

Zoom In Zoom Out Reset image size

For the purpose of assessing the resulting constraints on the DTD and λ we assume that BNS mergers from the current GW network will have precisely determined redshifts through associated EM counterparts and host galaxies. Therefore, the error budget is dominated by the Poisson error based on the detection rate. For the input model we assume Γ = −0.6, t_min = 700 Myr, and λ = 10⁻⁵ ${M}_{\odot }^{-1}$. Using our MCMC approach we show the resulting constraints on the DTD parameters for a year of AdLIGO/Virgo operations at design sensitivity in Figure 4. The results indicate that the DTD remains largely unconstrained, with the posterior distributions strongly influenced by the flat priors. In particular, t_min is unconstrained, while Γ and λ show a strong degeneracy, with median values that are biased away from the injected model. We find the same result for a decade of AdLIGO/Virgo operations.

Zoom In Zoom Out Reset image size

Our results for AdLIGO/Virgo suggest that even the proposed upgrades to the current facilities, such as A+ (Miller et al. 2015) and Voyager (Lantz et al. 2018) will not have a significant impact on the DTD as these facilities will still only detect BNS mergers in the local universe (see e.g., Figure 1, right panel of Reitze et al. 2019). As argued in Paper I, a more robust constraint on the DTD from the current generation of GW detectors may be achieved through the mass distribution of BNS merger host galaxies. However, even this approach leaves a lingering degeneracy between Γ and t_min.

The situation is drastically different for the third-generation detectors ET and CE. In Figure 5 we plot the expected detection rate as a function of redshift for ET, CE, and a network of ET+CE. Two improvements are readily apparent. First, the expected detection rate is about three orders of magnitude larger than for AdLIGO/Virgo. Second, the redshift range for BNS merger detections increases to z ∼ 5 in the case of ET+CE. The latter improvement results in a clear difference between the redshift distributions of the various DTDs, while the former improvement provides the detection statistics needed to distinguish between the DTD models. The differences between the various DTDs can no longer be scaled away with a change in λ (as is the case for AdLIGO/Virgo). In what follows we focus on the case of ET+CE as a realistic version of a third-generation detector network.

Zoom In Zoom Out Reset image size

In Figure 6 we show the result of MCMC fitting for the same input model used in the previous section. We show the results for 1 day, 1 week, 1 month, and 1 yr of observations. Unlike in the case of the current generation of GW detectors, we assume that the BNS mergers at cosmological distances will generally not have detectable EM counterparts (see Appendix C.). Instead we rely on distance information from the GW signal itself. We model the resulting redshift uncertainty as δz/z = 0.1z; a detailed motivation for this parameterization is provided in Appendix B.

Zoom In Zoom Out Reset image size

Our results show that an ET+CE network is able to constrain the DTD parameters and overcome the intrinsic degeneracy between Γ and t_min within a year of observations. The values of Γ and t_min can be determined to better than 10% accuracy. We note that these numbers depend on the overall event rate, which is determined by λ; here we use an injected value of 10⁻⁵ M_⊙⁻¹, but the results can be rescaled for higher or lower values.

To assess the impact of our input model on the results, in Figure 7 we repeat the same exercise, but for two different DTDs that favor short merger timescales: Γ = −1.2 with t_min = 30 Myr, and Γ = −1 with t_min = 100 Myr. Although the power-law index is recovered with the same accuracy as before, we find that t_min becomes more challenging to determine when its value is small. This is because the relative shift in the observed BNS merger redshift distribution becomes progressively smaller for small values of t_min, which is challenging to detect in the presence of realistic GW redshift uncertainties.

Zoom In Zoom Out Reset image size

We investigated how well the DTD and mass efficiency of BNS systems can be determined through the redshift distribution of BNS mergers detected by current and future GW networks. We model the DTD as a power law with a minimum merger timescale, and leave the mass efficiency as a free parameter. While other DTDs have been proposed (e.g., Simonetti et al. 2019), our primary conclusions should not be affected by the exact form of the DTD.

We find that current GW detectors, which can only detect BNS mergers in the local universe, cannot directly constrain the DTD due to their limited sensitivity. In effect, the various DTDs, with an appropriate scaling of λ, predict the same BNS merger detection rate at z ≈ 0. However, the situation is dramatically different for the anticipated third-generation detectors, which will be able to detect BNS mergers to z ≈ 5. For this cosmological population, even in the presence of redshift uncertainties of δz/z ≈ 0.1z from the GW data, the large detection rate and broad redshift range will precisely determine the DTD within about 1 yr of operations. It has been previously argued that the cosmological merger population uncovered by third-generation detectors will be able to constrain cosmological parameters (Sathyaprakash et al. 2010; Taylor et al. 2012; Taylor & Gair 2012; Vitale & Farr 2018); our results for the DTD further bolster the science case for third-generation GW detectors.

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. M.T.S. is thankful to Harvard-Smithsonian Center for Astrophysics for hospitality which made this work possible. S.V., K.K.Y.N., and C.W. acknowledge support of the National Science Foundation through the NSF award PHY1836814, the LIGO Laboratory, and the LIGO Data Grid clusters. LIGO was constructed by the California Institute of Technology and Massachusetts Institute of Technology with funding from the National Science Foundation and operates under cooperative agreement PHY-0757058.

The strain measured by a GW interferometer in frequency domain is given by

$$\begin{eqnarray}&&\tilde{h}(f)={F}_{+}{\tilde{h}}_{+}(f)+{F}_{\times }{\tilde{h}}_{\times }(f),\end{eqnarray} \tag{ 8 }$$

where ${\tilde{h}}_{+,\times }$ are the +, ×-polarization bases and F_+,× are the corresponding beam pattern functions,

$$\begin{eqnarray}&&{F}_{+}=g\left[\displaystyle \frac{1}{2}\left(1+{\cos }^{2}\theta \right)\cos 2\phi \cos 2\psi -\cos \theta \sin 2\phi \sin 2\psi \right],\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&{F}_{\times }=g\left[\displaystyle \frac{1}{2}\left(1+{\cos }^{2}\theta \right)\cos 2\phi \sin 2\psi +\cos \theta \sin 2\phi \cos 2\psi \right],\end{eqnarray} \tag{ 10 }$$

where θ, ϕ and ψ are the zenith, azimuth, and polarization angles, respectively, and g is a dimensionless coefficient determined by the geometry of an interferometer (Sathyaprakash & Schutz 2009; Schutz 2011; Chen et al. 2017). CE is a single interferometer with the angle between two arms equal to 90°, hence g_CE = 1 (Schutz 2011; Chen et al. 2017). ET consists of three identical interferometers with the angle between two arms equal to 60°, forming an equilateral triangle, hence ${g}_{\mathrm{ET}}=\sqrt{3}/2$ for each interferometer in ET (Punturo et al. 2010; Sathyaprakash et al. 2012).

The detection probability, P_det, is defined as the probability of a detection with ρ_net ≥ ρ_T, where ρ_T is the S/N threshold of detection and ${\rho }_{\mathrm{net}}^{2}={\sum }_{i}{\rho }_{i}^{2}$ is the network S/N as a geometric sum of S/N of each inteferometer. Assuming isotropic sky locations, orbital orientation, and polarization (i.e., uniform distribution of (cosθ, ϕ, ψ, cosι)), P_det is given analytically as

$$\begin{eqnarray}&&\begin{array}{l}{P}_{\det }\left({\theta }_{\mathrm{int}},z\right)\\ \quad \ =\ \displaystyle \int H\left(\displaystyle \frac{{\rho }_{T}}{{\rho }_{\mathrm{net}}({\theta }_{\mathrm{int}},z,\theta ,\phi ,\psi ,\iota )}\right)d\cos \theta d\phi d\psi d\cos \iota ,\end{array}\end{eqnarray} \tag{ 11 }$$

where θ_int is the set of intrinsic parameters, which are fixed at 1.4–1.4 M_⊙ and zero-spin for BNS systems, and H(w(z, θ, ϕ, ψ, ι)) is the unitary step function defined in w ∈ (0, 1] for w(z, θ, ϕ, ψ, ι) = ρ_T/ρ_net(z, θ, ϕ, ψ, ι).

For single CE or ET, we follow the inspiral approximation in Finn (1996). The S/N of a single interferometer is approximately

$$\begin{eqnarray}&&{\rho }^{2}\,=\,64{{\rm{\Theta }}}^{2}{\left(\displaystyle \frac{{r}_{0}}{{d}_{L}}\right)}^{2}{\left(\displaystyle \frac{{{ \mathcal M }}_{z}}{1.2{M}_{\odot }}\right)}^{5/3}{\zeta }^{2}({f}_{\max }),\end{eqnarray} \tag{ 12 }$$

where ${{\rm{\Theta }}}^{2}=4\left[{F}_{+}^{2}{\left(1+{\cos }^{2}\iota \right)}^{2}+4{F}_{\times }^{2}{\cos }^{2}\iota \right]$, and proportional to ${{\rm{\Theta }}}_{\mathrm{net}}^{2}={\sum }_{i}{{\rm{\Theta }}}_{i}^{2}$ for the same signal strain observed by homogeneous detectors such as single CE or ET. Hence P_det(w(z)) is equivalently the survival function of ${{\rm{\Theta }}}_{\mathrm{net}}/{{\rm{\Theta }}}_{\mathrm{net}}^{\max }$, where ${{\rm{\Theta }}}_{\mathrm{net}}^{\max }$ is the maximum response of a particular network of detectors (Finn 1996; Sathyaprakash & Schutz 2009; Schutz 2011; Chen et al. 2017). Here w(z) = ρ_T/ρ_opt(z), where ρ_opt is the optimal S/N at maximum Θ². CE has ${{\rm{\Theta }}}_{\mathrm{CE}}^{\max }=4$ and ET has ${{\rm{\Theta }}}_{\mathrm{ET}}^{\max }=4\times \sqrt{3\times 3/4}=6$. We have verified that the inspiral approximation in 3G detectors only results in few-percent difference in S/N for z ≲ 6, which is the region in which we are interested.

Approximated forms of P_det assuming uniformly distributed (cos θ, ϕ, ψ, cos ι) provided in Finn (1996) or Dominik et al. (2015) are only suitable for a single CE or a second-generation network. Therefore, we generate the survival function of Θ by drawing 10⁶ points of uniformly distributed (cos θ, ϕ, ψ, cos ι) in single CE or ET and fit the survival function with the following parametric form:

$$\begin{eqnarray}&&{P}_{\det }(w;A,B)=\displaystyle \frac{\mathrm{erf}\left(A-{Bw}\right)-\mathrm{erf}\left(A-B\right)}{\mathrm{erf}\left(A\right)-\mathrm{erf}\left(A-B\right)},\end{eqnarray} \tag{ 13 }$$

where $\mathrm{erf}(x)=\tfrac{2}{\sqrt{\pi }}{\int }_{0}^{x}{e}^{-{t}^{2}}{dt}$ is the error function. We also employ a tenth-order polynomial,

$$\begin{eqnarray}\begin{array}{rcl}{P}_{\det }(w;{a}_{k}) & = & \displaystyle \sum _{k=1}^{9}{a}_{k}{\left(1-w\right)}^{k}\\ & & +\ \left(1-\displaystyle \sum _{k=1}^{9}{a}_{k}\right){\left(1-w\right)}^{10},\end{array}\end{eqnarray} \tag{ 14 }$$

where {a_k} are the polynomial coefficients.

The above simplification using distribution of Θ breaks down for a CE+ET network due to the different sensitivity and heterogeneous geometry of each interferometer. Instead, we calculate the integral in Equation (11) by simulating waveform and network S/N for each redshift. Then w(z) = z/z_horizon and z_horizon is the redshift of the detector horizon, which is ∼12.5 in CE+ET for a BNS merger. Again we fit the simulation result for CE+ET using Equations (13) and 14.

Tables 1 and 2 show the definitions of w and fitting parameters of Equations (13) and 14 in each network. Figure 8 shows the comparison of actual simulation and the parametric fits.

Zoom In Zoom Out Reset image size

In the case of third-generation detectors it is unlikely that most BNS merger detections will have EM counterparts. Instead, the distance information will need to be gleaned from the GW signal itself. There is an inherent degeneracy between a binary's inclination angle in the sky with respect to a GW detector, θ_JN, and the luminosity distance (Schutz 2011; Abbott et al. 2016; Chen et al. 2018; Usman et al. 2018). Messenger & Read (2012) show that for a range of representative neutron star equations of state the redshift of such systems can be determined to an accuracy of ∼8%–40% for z < 1 and ∼9%–65% for 1 < z < 4. However, for a binary to be detectable at high redshifts, we expect the inclination to be close to face-on (θ_JN = 0°), or face-off (θ_JN = 180°), as most of the energy in GWs is released along the angular momentum vector of the binary. Schutz (2011) derived an analytic formulation for the distribution of inclination angle of sources detectable by advanced detectors; ≲7% (≲3%) of detectable events will have viewing angles of >70° (>80°) (Chen et al. 2018).

In principle, we can simulate BNS signals at different redshifts in third-generation detectors and estimate the evolution of uncertainty in parameter estimation using MCMC. Because BNSs are low-mass systems with a long coalescence time, the parameter estimation is computationally expensive. To approximate the redshift uncertainty of a BNS at fixed redshift, we may extrapolate a BBH signal by lowering the starting frequency from 10 to 5 Hz, as to mimic the long inspiral phase in a BNS merger. We obtain the redshift uncertainties of BBHs in CE+ET from the simulations of Vitale & Whittle (2018), and fit the mean redshift uncertainties as a function of true redshifts. Then we calibrate our fit using the above extrapolation scheme to approximate the redshift uncertainties as δz/z ≈ 0.1z.

Alternatively, Chen et al. (2018) developed a rapid algorithm that provides a luminosity distance uncertainty estimate for a large population of BNS merger detections. We use this algorithm to simulate 2000 BNS detections in a third-generation network, and compare the results to the extrapolation procedure above. We find that the two approaches yield consistent results. We therefore use the extrapolated distance uncertainty, and convert it to the redshift uncertainty through the adopted cosmology in this work.

The distance-inclination degeneracy can be broken through the detections of an associated SGRB. The relative fraction of on-axis mergers is only a few percent. A γ-ray detection alone will thereby reduce the overall redshift uncertainty of at most a few percent of BNS merger detections, and likely over a restricted redshift range (perhaps to z ∼ 2). An afterglow detection can further lead to a precise redshift determination through an associated host galaxy, but to date such detections have mainly been limited to z ≲ 1 (Berger 2014), which is generally not high enough to make an impact on the DTD determination (see Figure 5). Thus, it seems unlikely that associated SGRBs and their afterglows will substantially improve the constraints on the DTD.