Future Prospects for Ground-based Gravitational-wave Detectors: The Galactic Double Neutron Star Merger Rate Revisited

Since we want to calculate the total number of DNS systems like those that have been observed, we fix the physical parameters of the pulsars generated in our simulation to represent the DNS systems in which we are interested. These physical parameters include the pulse period, pulse width, and orbital parameters like eccentricity, orbital period, and semimajor axis.

However, even if the physical parameters of the pulsars are the same, their luminosity will not be the same. Thus, to model the luminosity distribution of these pulsars, we use a log-normal distribution with a mean of $\left\langle {\mathrm{log}}_{10}L\right\rangle =-1.1$ (L = 0.07 mJy kpc²) and standard deviation, ${\sigma }_{{\mathrm{log}}_{10}L}=0.9$ (Faucher-Giguère & Kaspi 2006).

We also vary the spectral index of the simulated pulsar population. We assume the spectral indices have a normal distribution, with mean α = −1.4, and standard deviation β = 1 (Bates et al. 2013).

All of the DNS systems that merge within a Hubble time have either been detected or discovered in the following surveys: the Pulsar Arecibo L-band Feed Array survey (PALFA; Cordes et al. 2006), the High Time-Resolution Universe pulsar survey (HTRU; Keith et al. 2010), the Parkes High-latitude pulsar survey (Burgay et al. 2006), the Parkes Multibeam Survey (Manchester et al. 2001), and the survey carried out by Wolszczan (1991) in which B1534+12 was discovered. All of these surveys together cover more area on the sky than that covered by the 18 surveys simulated in KKL03 and by Kim et al. (2015), who added the Parkes Multibeam Pulsar survey to the latter.

We implement these surveys in our simulations with PsrPopPy. We generate a survey file (see Section 4.1 in Bates et al. 2014) for each of them using the published survey parameters. These parameters are then used to estimate the radiometer noise in each survey which, along with a fiducial signal-to-noise cut-off, will determine whether a pulsar from the simulated population can be detected with a given survey. For example, one important difference in these surveys is their integration time, which ranges from 34 s for Arecibo drift-scan surveys to 2100 s in the Parkes Multibeam survey. Other selection effects can be introduced through differences in the sensitivity of the different surveys, the portion and area of the sky covered, and minimum signal-to-noise ratio cut-offs.

PSR J1757–1854 was discovered in the HTRU low-latitude survey using a novel search technique (Cameron et al. 2018). As described in Ng et al. (2015), the original integration time of 4300 s of the HTRU survey was successively segmented by a factor of two into smaller time intervals until a pulsar was detected. This had the effect of reducing Doppler smearing due to extreme orbital motion in tight binary systems (see Section 2.4 for more on Doppler smearing). The shortest segmented integration time used in their analysis was 537 s (one-eighth segment), which implies that the data are sensitive to binary systems with orbital periods ${P}_{{\rm{b}}}\geqslant 1.5\,\mathrm{hr}$ (Ng et al. 2015). All of these segments are searched for pulsars in parallel. We use the integration time of 537 s in our analysis to ensure that the HTRU survey is sensitive to all the DNS systems included in this analysis. We demonstrate the effect of this choice in Section 4.2.1.

The survey files are available in the GitHub repository associated with this paper.

For the radial distribution of the DNS systems in the Galaxy, like Kim et al. (2015) we use the model proposed in Lorimer et al. (2006). For the distribution of pulsars in terms of their height, z, with respect to the Galactic plane, we use the standard two-sided exponential function (Lyne 1998; Lorimer et al. 2006),

$$\begin{eqnarray}&&f(z)\propto \exp \left(\displaystyle \frac{-| z| }{{z}_{0}}\right)\end{eqnarray} \tag{ 2 }$$

where z₀ is the vertical scale height. To constrain z₀, we simulate DNS populations with a uniform period distribution ranging from 15 to 70 ms, consistent with the periods of the recycled pulsars in the DNS systems listed in Table 1, and the aforementioned luminosity and spectral index distribution. We generate these populations with vertical scale heights ranging from ${z}_{0}=0.1\,\mathrm{kpc}$ to ${z}_{0}=2\,\mathrm{kpc}$. We run the surveys described in Section 2.2 on these populations to determine the median vertical scale height of the pulsars that are detected in these surveys. We also calculate the median $\mathrm{DM}\times \sin (| b| )$, which is more robust against errors in converting from dispersion measure to a distance using the NE2001 Galactic electron density model (Cordes & Lazio 2002).

We compare these values at different input vertical scale heights with the corresponding median values for the real DNS systems. We show the median $\mathrm{DM}\times \sin (| b| )$ value and the median vertical z-height of the pulsars detected in the simulations as a function of the input z₀ in Figures 1 and 2 respectively. In both of these plots, the median values of the real DNS population are plotted as the red dashed line, with the error on the median shown by the shaded cyan region. As can be seen, the analysis using $\mathrm{DM}\times \sin (| b| )$ predicts a vertical scale height of z₀ = 0.4 ± 0.1 kpc, while the analysis using the z-height estimated using the NE2001 model (Cordes & Lazio 2002) returns a vertical scale height of ${z}_{0}={0.4}_{-0.2}^{+0.3}$ kpc. While both these values are consistent with each other, the vertical scale height returned by the $\mathrm{DM}\times \sin (| b| )$ analysis yields a better constraint on the scale height which is more in line with vertical scale heights for ordinary pulsars (0.33 ± 0.03 kpc; Mdzinarishvili & Melikidze 2004; Lorimer et al. 2006) and millisecond pulsars (0.5 kpc; Levin et al. 2013). We expect the neutron stars that exist in DNS systems, and particularly those DNS systems that merge within a Hubble time, to be born with small natal kicks so as not to disrupt the orbital system. Consequently, we would expect these systems to be closer to the Galactic plane than the general millisecond pulsar population. As a result, we adopt a vertical scale height of z = 0.33 kpc as a conservative estimate on the vertical scale height of the DNS population distribution. This difference in the vertical scale height does not result in a significant change in the merger rates.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The motion of the pulsar in the orbit of the DNS system introduces a Doppler shifting of the observed pulse period. The extent of the Doppler shift depends on the velocity and acceleration of the pulsar in different parts of its orbit. This Doppler shift results in a reduction in the signal-to-noise ratio for the observation of the pulsar (Bagchi et al. 2013).

To account for this effect, we use the algorithm developed by Bagchi et al. (2013), which quantifies the reduction in the signal-to-noise ratio as a degradation factor, $0\lt \gamma \lt 1$, averaged over the entire orbit. This degradation factor depends on the orbital parameters of the DNS system (such as eccentricity and orbital period), the mass of the two neutron stars, the integration time for the observation, and the search technique used in the survey (for example, HTRU and PALFA use acceleration-searches; Bagchi et al. 2013). A degradation factor $\gamma \sim 1$ implies very little Doppler smearing, while a degradation factor $\gamma \sim 0$ implies heavy Doppler smearing in the pulsar's radio emission.

The implementation of the algorithm as a Fortran program was kindly provided to us by the authors of Bagchi et al. (2013), which we make available⁵ with their permission. Since PsrPopPy does not include functionality to handle this degradation factor, we had to manually introduce it into the source code of PsrPopPy. The modified PsrPopPy source files are also available on the GitHub repository.

The beaming correction factor, f_b, is defined as the inverse of the pulsar's beaming fraction, i.e., the solid angle swept out by the pulsar's radio beam divided by $4\pi$. PSRs B1913+16, B1534+12, and J0737–3039A/B have detailed polarimetric observation data, from which precise measurement of their beaming fractions, and thus their beaming correction factors, has been possible. These beaming corrections are collected in Table 2 of Kim et al. (2015).

However, the other merging DNS systems are relatively new discoveries and do not have measured values for their beaming fractions. Thus, we assume that the beaming correction factor for these new pulsars is the average of the measured beaming correction factor for the three aforementioned pulsars, i.e., 4.6. We list these beaming fractions in Table 2, and defer discussion of their effect on the merger rate for Section 4.

Table 2.  Parameters and Results from the DNS Merger Rate Analysis Described in Section 3

PSR	fb	δ	τage	Nobs	Npop	${ \mathcal R }$
			(Myr)			(Myr−1)
B1534+12	6.0	0.04	208	${{\bf{114}}}_{-{\bf{80}}}^{+{\bf{519}}}$	${{\bf{688}}}_{-{\bf{483}}}^{+{\bf{2546}}}$	${\bf{0}}.{{\bf{2}}}_{-{\bf{0.1}}}^{+{\bf{1.1}}}$
J1756−2251	4.6	0.03	396	${{\bf{138}}}_{-{\bf{96}}}^{+{\bf{627}}}$	${{\bf{633}}}_{-{\bf{440}}}^{+{\bf{2878}}}$	${\bf{0}}.{{\bf{4}}}_{-{\bf{0.3}}}^{+{\bf{1.7}}}$
J1913+1102	4.6	0.06	2625	${{\bf{182}}}_{-{\bf{132}}}^{+{\bf{830}}}$	${{\bf{834}}}_{-{\bf{605}}}^{+{\bf{3813}}}$	${\bf{0}}.{{\bf{3}}}_{-{\bf{0.2}}}^{+{\bf{1.2}}}$
J1906+0746	4.6	0.01	0.11	${{\bf{62}}}_{-{\bf{36}}}^{+{\bf{276}}}$	${{\bf{284}}}_{-{\bf{165}}}^{+{\bf{1265}}}$	${\bf{4}}.{{\bf{7}}}_{-{\bf{2.7}}}^{+{\bf{21.1}}}$
B1913+16	5.7	0.169	77	${{\bf{182}}}_{-{\bf{132}}}^{+{\bf{822}}}$	${{\bf{1040}}}_{-{\bf{754}}}^{+{\bf{4706}}}$	${\bf{2}}.{{\bf{8}}}_{-{\bf{2.0}}}^{{\bf{12.3}}}$
J0737−3039A	2.0	0.27	159	${{\bf{465}}}_{-{\bf{347}}}^{+{\bf{2097}}}$	${{\bf{931}}}_{-{\bf{695}}}^{+{\bf{4193}}}$	${\bf{3}}.{{\bf{9}}}_{-{\bf{2.9}}}^{+{\bf{17.4}}}$
J1757−1854	4.6	0.06	87	${{\bf{198}}}_{-{\bf{144}}}^{+{\bf{906}}}$	${{\bf{908}}}_{-{\bf{660}}}^{+{\bf{4161}}}$	${\bf{5}}.{{\bf{6}}}_{-{\bf{4.1}}}^{+{\bf{25.7}}}$
J1946+2052	4.6	0.06	247	${{\bf{306}}}_{-{\bf{224}}}^{+{\bf{1397}}}$	${{\bf{1402}}}_{-{\bf{1026}}}^{+{\bf{6417}}}$	${\bf{4}}.{{\bf{8}}}_{-{\bf{5.5}}}^{+{\bf{21.9}}}$

Note. Here, f_b is the beaming correction factor, δ is the pulse duty cycle, τ_age is the effective age described in Section 3, N_obs is the number of each DNS system that are beaming toward the Earth, N_pop is the total number of each DNS system in the Milky Way, and ${ \mathcal R }$ is the merger rate of each individual DNS system population, the probability distribution function for which is shown in Figure 3. Errors are quoted at the 95% confidence interval.

Download table as:  ASCIITypeset image

The effective lifetime of a DNS binary, ${\tau }_{\mathrm{life}}$, is defined as the time interval during which the DNS system is detectable. Thus, it is the sum of the time since the formation of the DNS system and its remaining lifetime,

$$\begin{eqnarray}\begin{array}{rcl}{\tau }_{\mathrm{life}} & = & {\tau }_{\mathrm{age}}+{\tau }_{\mathrm{obs}}\\ & = & \min \left({\tau }_{{\rm{c}}},{\tau }_{{\rm{c}}}{\left[1-\displaystyle \frac{{P}_{\mathrm{birth}}}{{P}_{{\rm{s}}}}\right]}^{n-1}\right)+\min ({\tau }_{\mathrm{merger}},{\tau }_{{\rm{d}}}).\end{array}\end{eqnarray} \tag{ 3 }$$

Here ${\tau }_{{\rm{c}}}={P}_{{\rm{s}}}/(n-1)\dot{{P}_{{\rm{s}}}}$ is the characteristic age of the pulsar, n is the braking index, assumed to be 3, ${P}_{\mathrm{birth}}$ is the period of the millisecond pulsar at birth, i.e., when it begins to move away from the fiducial spin-up line on the $P-\dot{P}$ diagram, ${P}_{{\rm{s}}}$ is the current spin period of the pulsar, ${\tau }_{\mathrm{merger}}$ is the time for the DNS system to merge, and ${\tau }_{{\rm{d}}}$ is the time in which the pulsar crosses the "death line" beyond which it should not radiate significantly (Chen & Ruderman 1993).

Unlike normal pulsars, the characteristic age, ${\tau }_{{\rm{c}}}$, for millisecond and recycled pulsars may not be a very good indicator of their true age. This is due to the fact that the period of the pulsar at birth is much smaller than its current period, which is not true for recycled millisecond pulsars found in DNS systems. A better estimate for the age of a recycled millisecond pulsar can be calculated by measuring its distance from a fiducial spin-up line on the $P-\dot{P}$ diagram (Arzoumanian et al. 1999), represented by the second part of the first term in Equation (3).

Finally, the time for which a given DNS system is detectable after birth depends on whether we are observing the non-recycled companion pulsar (J0737–3039B, J1906+0746) or the recycled pulsar in the DNS system (e.g., B1913+16, J1757–1854, J1946+2052, etc.). In the latter case, the combination of a small spin-down rate and millisecond period ensures that the DNS system remains detectable until the epoch of the merger. However, for the former case, both the period and spin-down rate are at least an order of magnitude larger than their recycled counterparts. As such, the time for which these systems are detectable depends on whether they cross the pulsar "death line" before their epoch of merger (Chen & Ruderman 1993). The radio lifetime of any pulsar is defined as the time it takes the pulsar to cross this fiducial "death line" on the $P-\dot{P}$ diagram (Chen & Ruderman 1993).

We estimate the radio lifetime for J1906+0746 using two different techniques. One estimate is described by Chen & Ruderman (1993) and assumes a simple dipolar rotator to find the time to cross the death line. Using Equation (6) in their paper we calculate a radio lifetime of ${\tau }_{{\rm{d}}}\sim 3\,\mathrm{Myr}$ for J1906+0746. However, as discussed in Chen & Ruderman (1993), the death line for a pure dipolar rotator might not be an accurate turn-off point for pulsars, with many observed pulsars lying past this line on the $P-\dot{P}$ diagram. A better estimate of the radio lifetime might be given by Equation (9) in Chen & Ruderman (1993), which assumes a twisted field configuration for pulsars. Using this, we find ${\tau }_{{\rm{d}}}\sim 30\,\mathrm{Myr}$. Another estimate for the radio lifetime can be made from spin-down energy loss considerations. Adopting the formalism given in van den Heuvel & Lorimer (1996), we find, for a simple dipolar spin-down model, that a pulsar with a current spin-down energy loss rate ${\dot{E}}_{\mathrm{now}}$ and characteristic age τ will reach a cut-off $\dot{E}$ value of 10³⁰ erg s⁻¹ below which radio emission through pair production is suppressed on a timescale:

$$\begin{eqnarray}&&{\tau }_{{\rm{d}}}={\tau }_{{\rm{c}}}\left(\sqrt{\displaystyle \frac{{\dot{E}}_{\mathrm{now}}}{\dot{E}={10}^{30}}}-1\right).\end{eqnarray} \tag{ 4 }$$

Using this formalism, we calculate a radio lifetime of ${\tau }_{{\rm{d}}}\sim 60\,\mathrm{Myr}$. This method has been used in previous estimates (Kalogera et al. 2001; Kim et al. 2003, 2015) of the merger rates, and represents a conservative estimate on the radio lifetime of J1906+0746. We adopt it here as the fiducial radio lifetime of J1906+0746 for consistency, and defer the discussion of the implications of variation in the calculated radio lifetime to Section 4.

A similar analysis could be done for pulsar B in the J0737–3039 system. However, unlike Kim et al. (2015), we do not include B in our merger rate calculations. The uncertainties in the radio lifetime are very large, as for PSR J1906+0746, and therefore pulsar A provides a much more reliable estimate of the numbers of such systems. In addition, unlike J1906+0746, pulsar B also shows large variations in its equivalent pulse width (Kim et al. 2015), and thus its duty cycle, due to pulse profile evolution through geodetic precession (Perera et al. 2010). This also leads to an uncertainty in its beaming correction factor (see Figure 4 in Kim et al. 2015). There are additional uncertainties introduced by pulsar B exhibiting strong flux density variations over a single orbit around A. All these factors introduce a large uncertainty in the merger rate contribution from B, and do not provide better constraints on the merger rate compared to when only pulsar A is included (Kim et al. 2015). Finally, the Double Pulsar system was discovered through pulsar A and will remain detectable through pulsar A long after B crosses the death line. For these reasons, we do not include pulsar B in our analysis.

After calculating individual merger rates from each DNS system, we need to combine these merger rate probability distributions to find the combined Galactic probability distribution. We can do this by treating the merger rate for the individual DNS systems as independent continuous random variables. In that case, the total merger rate for the Galaxy will be the arithemtic sum of the individual merger rates

$$\begin{eqnarray}&&{{ \mathcal R }}_{\mathrm{MW}}=\displaystyle \sum _{i=1}^{8}{{ \mathcal R }}_{i}\end{eqnarray} \tag{ 12 }$$

with the total Galactic merger rate probability distribution given by a convolution of the individual merger rate probability distributions,

$$\begin{eqnarray}&&P({{ \mathcal R }}_{\mathrm{MW}})=\displaystyle \prod _{i=1}^{8}P({{ \mathcal R }}_{i})\end{eqnarray} \tag{ 13 }$$

where $\prod$ denotes convolution. As the number of known DNS systems increases over time, the method of convolution of individual merger rate probability density functions is more efficient than computing an explicit analytic expression as in KKL03 and Kim et al. (2015).

Combining all the individual Galactic merger rates, we obtain a total Galactic merger rate of ${{ \mathcal R }}_{\mathrm{MW}}={42}_{-14}^{+30}$ Myr⁻¹, which is shown in Figure 4.

Zoom In Zoom Out Reset image size

The Galactic merger rate calculated above can be extrapolated to calculate the number of DNS merger events that LIGO will be able to detect. Assuming that the DNS formation rate is proportional to the formation rate of massive stars, which is in turn proportional to the B-band luminosity of a given galaxy (Phinney 1991; Kalogera et al. 2001), the DNS merger rate within a sphere of radius D is given by (Kopparapu et al. 2008)

$$\begin{eqnarray}&&{{ \mathcal R }}_{\mathrm{LIGO}}={{ \mathcal R }}_{\mathrm{MW}}\left(\displaystyle \frac{{L}_{\mathrm{total}}(D)}{{L}_{\mathrm{MW}}}\right)\end{eqnarray} \tag{ 14 }$$

where L_total(D) is the total blue luminosity within a distance D, and ${L}_{\mathrm{MW}}=1.7\times {10}^{10}{L}_{B,\odot }$, where ${L}_{B,\odot }=2.16\times {10}^{33}$ erg s⁻¹, is the B-band luminosity of the Milky Way (Kopparapu et al. 2008).

Using a reference LIGO range distance of ${D}_{{\rm{r}}}=100$ Mpc (Abbott et al. 2018), and following the arguments laid out in Kopparapu et al. (2008), we can calculate the rate of DNS merger events visible to LIGO (Equation (19) in Kopparapu et al. 2008):

$$\begin{eqnarray}\begin{array}{rcl}{{ \mathcal R }}_{\mathrm{LIGO}} & = & \displaystyle \frac{N}{T}\\ & = & 7.4\,\times \,{10}^{-3}\left(\displaystyle \frac{{ \mathcal R }}{{\left({10}^{10}{L}_{B,\odot }\right)}^{-1}\,{\mathrm{Myr}}^{-1}}\right)\\ & & \times \,{\left(\displaystyle \frac{{D}_{{\rm{r}}}}{100\mathrm{Mpc}}\right)}^{3}{\mathrm{yr}}^{-1},\end{array}\end{eqnarray} \tag{ 15 }$$

where N is the number of mergers in T years, and ${ \mathcal R }={{ \mathcal R }}_{\mathrm{MW}}/{L}_{\mathrm{MW}}$ is the Milky Way merger rate weighted by the Milky Way B-band luminosity.

Using the above equation, we calculate the DNS merger detection rate for LIGO,

$$\begin{eqnarray}&&{{ \mathcal R }}_{\mathrm{LIGO}}\equiv {{ \mathcal R }}_{\mathrm{PML}18}={0.18}_{-0.06}^{+0.13}\times {\left(\displaystyle \frac{{D}_{{\rm{r}}}}{100\mathrm{Mpc}}\right)}^{3}{\mathrm{yr}}^{-1},\end{eqnarray} \tag{ 16 }$$

where we use ${{ \mathcal R }}_{\mathrm{PML}18}$ to distinguish our merger detection rate estimate from the others that will be referred to later in the paper.

The recent detection of a DNS merger by LIGO (Abbott et al. 2017a) enabled a calculation of the rate of DNS mergers visible to LIGO. The rate that was calculated by Abbott et al., converted to the units used in our calculations, is

$$\begin{eqnarray}&&{{ \mathcal R }}_{\mathrm{LIGO}}\equiv {{ \mathcal R }}_{{\rm{A}}17}={1.54}_{-1.22}^{+3.20}\times {\left(\displaystyle \frac{{D}_{{\rm{r}}}}{100\mathrm{Mpc}}\right)}^{3}{\mathrm{yr}}^{-1}\end{eqnarray} \tag{ 17 }$$

where ${{ \mathcal R }}_{{\rm{A}}17}$ is the merger detection rate and the errors quoted are 90% confidence intervals.

We plot both the rate estimates in Figure 5. This rate estimated by LIGO is in agreement with the DNS merger detection rate that we calculate using the Milky Way DNS binary population, ${{ \mathcal R }}_{\mathrm{PML}18}$, at the upper end of the 90% confidence level range.

Zoom In Zoom Out Reset image size

4.2.1. Luminosity Distribution

4.2.2. Beaming Correction Factors

4.2.3. The Effective Lifetime of J1906+0746

4.2.4. Extrapolation to LIGO's Observable Volume

4.2.5. Unobserved Underlying DNS Population in the Milky Way

We can also compare our merger detection rate to that predicted through theoretical studies and simulations of the formation and evolution of DNS binary systems. This approach to calculating the merger detection rate factors in the different evolutionary scenarios leading to the formation of the DNS system, including modeling stellar wind in progenitor massive star binaries, core-collapse and electron-capture supernovae explosions, natal kicks to the NSs, and the common-envelope phase (Abadie et al. 2010; Dominik et al. 2013, 2015). We compare our merger detection rate to the predictions made using the above methodology following the DNS merger detected by LIGO (Abbott et al. 2017a), i.e., the studies by Chruslinska et al. (2018) and Mapelli & Giacobbo (2018). We plot their estimates along with those calculated in this work in Figure 6.

Zoom In Zoom Out Reset image size

Chruslinska et al. (2018), using their reference model, calculated a merger detection rate density for LIGO of 48.4 Gpc⁻³ yr⁻¹ which, scaled to a range distance of 100 Mpc, is equivalent to a merger detection rate of 0.0484 yr⁻¹. This value is significantly lower than our range of predicted merger detection rates. In addition to the reference model, they also calculate the merger detection rate densities for a variety of different models, with the most optimistic model predicting a value of ${600}_{-300}^{+600}$ Gpc⁻³ yr⁻¹. Scaling this to our reference range distance of 100 Mpc, we obtain a merger detection rate of ${0.6}_{-0.3}^{+0.6}$ yr⁻¹, which is consistent with the LIGO calculated value (${{ \mathcal R }}_{{\rm{A}}17}={1.54}_{-1.22}^{+3.20}$ yr⁻¹). However, this optimistic model assumes that Hertzsprung gap donors avoid merging with their binary companions during the common-envelope phase. Applying the same evolutionary scenario to black hole binaries (BHBs) overestimates their merger detection rate from that derived using the BHB mergers observed by LIGO (Chruslinska et al. 2018). Thus, for the optimistic model to be correct would need the common-envelope process to work differently for BHB systems as compared to DNS systems, or that BHB systems would endure a bigger natal kick in the same formation scenario than DNS systems would (Chruslinska et al. 2018).

Mapelli & Giacobbo (2018) showed that the above problem could be avoided and a rate consistent with the LIGO prediction of the merger detection rate (${{ \mathcal R }}_{{\rm{A}}17}={1.54}_{-1.22}^{+3.20}$ yr⁻¹) could be obtained if there is high efficiency in the energy transfer during the common-envelope phase coupled with low kicks for both electron capture and core-collapse supernovae. Based on their population synthesis, they calculate a merger detection rate density of ∼600 Gpc⁻³ yr⁻¹ (for α = 5, low σ in Figure 1 of Mapelli & Giacobbo 2018). The full range of merger detection rate densities predicted by Mapelli & Giacobbo ranges from ∼20 Gpc⁻³ yr⁻¹ to ∼600 Gpc⁻³ yr⁻¹, which at a range distance of 100 Mpc corresponds to a merger detection rate ranging from 0.02 yr⁻¹ to 0.6 yr⁻¹. This rate is consistent with that derived in this work, and lends credence to the hypotheses of high energy transfer efficiency in the common-envelope phase and low natal kicks in DNS systems made by Mapelli & Giacobbo (2018).

In the short term, the difference between our merger rate and that calculated by LIGO can be clarified from the results of the third operating run (O3), which is scheduled for early 2019. Based on the fiducial model in our analysis and the predicted range distance of 120–170 Mpc for O3 (Abbott et al. 2018), we predict, accounting for 90% confidence intervals, that LIGO–Virgo will detect anywhere between zero and two DNS mergers. Further detections or non-detections by LIGO will be able to shed light on the detection rate within LIGO's observable volume. In addition, the localization of these mergers to their host galaxies as demonstrated by GW170817 (Abbott et al. 2017b) will determine the contribution of galaxies lacking in blue luminosity (such as ellipticals) to the total merger rate.

In the long term, with the advent of new large-scale telescope facilities such as the Square Kilometer Array (Carilli & Rawlings 2004), we should be able to survey our Galaxy with a much higher sensitivity. Such deep surveys might reveal more of the DNS population, which would yield a better constraint on the Galactic merger rate.

In addition to future radio surveys, a large number of LIGO detections of DNS mergers will allow us to probe the underlying DNS population directly. Assuming no large deviations from the DNS population parameters adopted in this study (see Sections 2 and 4), a significantly larger number of DNS merger detections by LIGO would imply a larger underlying DNS population. The localization of the DNS mergers to their host galaxies will allow us to test the variation in the DNS population with respect to host galaxy morphology. We might also be able to test if the DNS population in different galaxies is similar to that in the Milky Way. This will clarify the effect of the host galaxy morphology on the evolutionary scenario of DNS systems.