SHATTERING FLARES DURING CLOSE ENCOUNTERS OF NEUTRON STARS

Kocsis et al. (2006) calculated the parabolic encounter rate for compact objects in globular clusters using simplified globular cluster models, predicting a rate of ≳ 1 detection per year for advanced LIGO in optimistic scenarios. However, their detection rates are dominated by rare distant events involving close encounters of ≳ 20 M_☉ BHs.

Lee et al. (2010) examine various dynamical pathways to sGRBs in globular clusters, including binary interactions and tidal capture. They calculate a high rate of close encounters for two NSs, $\Gamma _{\rm NS-NS}^{({\rm GC})} \sim 55$ yr⁻¹ Gpc⁻³, using as a calibration the estimate of ∼10⁴ NSs in the collapsed core of M15, from the Fokker–Planck calculations of Dull et al. (1997). This would require an extremely high NS retention fraction. Subsequent more careful calculations by Murphy et al. (2011) have determined the number of NSs in the core of M15 to be closer to ∼10³, which is consistent with ∼1%–10% retention fraction estimates from pulsar kick velocities (Drukier 1996; Hansen & Phinney 1997; Davies & Hansen 1998). This reduces the estimates of Lee et al. (2010) to $\Gamma _{\rm NS-NS}^{({\rm GC})} \sim 0.5$ yr⁻¹ Gpc⁻³, however, Samsing et al. (2013) have recently showed that considering chaotic resonances in binary–single interactions can significantly increase the expected rate of eccentric binaries.

While there are many globular clusters per galaxy, high kick velocities at NS birth significantly lower their retention fraction. The deeper potentials of galactic nuclei may provide dense stellar environments where the NS retention fraction is higher and close encounters are more likely to occur. O'Leary et al. (2009) and Kocsis & Levin (2012) both provide estimates of ∼1–100 BH close passages per year within a few Gpc detectable by advanced LIGO, with GW-detectable NS– BH encounters estimated to be ∼1% of this. However, as we discuss in Appendix A below, they scale by a factor ξ, representing the contribution due to the variance of the nuclear cluster density across galaxies. They take this factor to be ξ ≳ 30–100, but we find that ξ is more correctly evaluated to be significantly lower, even with the most optimistic assumptions.

In Appendix A, we have re-evaluated the rates for single–single eccentric captures of compact objects in nuclear star clusters containing massive central BHs, assuming a simplified isothermal density distribution, as in Kocsis & Levin (2012). We find that the 10 M_☉ BH–BH eccentric capture encounter rate is $\Gamma _{\rm BH-BH}^{({\rm GN}, {\rm EC})}\sim 0.02$ yr⁻¹ Gpc⁻³, which is significantly lower than the previously estimated values. However, a more top-heavy initial mass function (IMF; Bartko et al. 2010), along with enhanced segregation and spatial flattening for heavier BHs, may help to increase this value.

The NS–NS and NS–BH (10 M_☉) rates can be estimated in a similar fashion for an isothermal ∝r⁻² density distribution to be $\Gamma _{\rm NS-NS}^{({\rm GN},{\rm EC})} \sim 0.04\hbox{--}6$ yr⁻¹ Gpc⁻³, and $\Gamma _{\rm NS-BH}^{({\rm GN},{\rm EC})} \sim 0.05\hbox{--} 0.6$ yr⁻¹ Gpc⁻³ with the range mainly due to uncertainty in the IMF and mass loss models for NS progenitors (O'Connor & Ott 2011). However, there is reason to suspect that the slope of the density distribution is somewhat flattened due to interaction with segregated BHs (see, e.g., O'Leary et al. 2009 and references therein). Taking a NS density distribution ∝r^−3/2 as a lower bound for our rates, we find in this case $\Gamma ^{({\rm GN},{\rm EC})}_{\rm NS-NS} \simeq 0.003\hbox{--} 0.3$ yr⁻¹ Gpc⁻³.

Note that in evaluating the above rates, we have made very optimistic assumptions about systematic versus intrinsic variation of the relevant observations, as in O'Leary et al. (2009). Possibly more realistic estimates for this intrinsic scatter reduce these rates by a factor of ∼4.

In the high density, high relative-velocity region near the center of nuclear star clusters, there may be a significant rate of hyperbolic passages where the periapse is sufficiently close to trigger a shattering flare, but insufficient to result in eccentric capture by GW emission. In Appendix B, we have calculated the rate of encounters that result in a shattering flare during the first close passage, and find this to be higher than the eccentric capture rate for the fiducial model used. We find, for the optimistic assumptions about intrinsic variation used above, and assuming our canonical isothermal model, $\Gamma ^{({\rm GN}, {\rm SF})}_{\rm NS-NS} \simeq 0.2\hbox{--} 60$ yr⁻¹ Gpc⁻³. For a more flattened density profile ∝r^−3/2, we have $\Gamma ^{({\rm GN}, {\rm SF})}_{\rm NS-NS} \simeq 0.005\hbox{--}0.5$ yr⁻¹. Taking our less generous estimates for the intrinsic scatter across galaxies significantly reduces these rates by a factor of ∼6.

In the above discussion, we have primarily considered single–single interactions in determining the event rates in dense clusters. Binary–single or binary–binary interactions have larger cross sections and could increase the rates for such events (see, e.g., Miller & Lauburg 2009). Recently, Katz & Dong (2012) and Kushnir et al. (2013) demonstrated that Kozai–Lidov type interactions can drive the inner binaries of hierarchical triple systems toward extreme eccentricity, with collisions occurring when the periapse distance is driven below the stellar radius. They claim that the rates for such Kozai-oscillation driven collisions between white dwarfs in field triple star systems can be comparable to the SN Type Ia rate. Similar interactions could potentially drive close encounters of NSs or BHs in triple systems. However, the periapse for shattering flares is two orders of magnitude smaller than those considered by Katz & Dong (2012), and NSs and BHs are substantially more rare than white dwarfs, particularly outside of dense clusters. Within clusters, such triple systems would need to reach these extreme eccentricities quickly, before other encounters ionize away the softer less-bound outer companion.

The cross section for eccentric capture is basically the cross section for which the energy emitted by GW (or lost due to tidal interactions) during a close encounter exceeds the kinetic energy of the objects at infinity. The maximum periapse for capture is given by Quinlan & Shapiro (1989) as

$$\begin{eqnarray} r_{p, {\rm max}} & = & \left[\frac{85 \pi \sqrt{2} G^{7/2} m_i m_j (m_i + m_j)^{3/2}}{12 c^5 |{\boldsymbol v}_i - {\boldsymbol v}_j|^2} \right]^{2/7} \nonumber\\ &=& 190\, {\rm km} \left(\frac{\eta }{0.25} \right)^{2/7} \left(\frac{m_{{\rm tot}}}{2.8 \ M_\odot } \right) \left(\frac{v_{\rm rel}}{10^3\ {\rm km\, s}^{-1}}\right)^{-4/7},\nonumber\\ \end{eqnarray} \tag{ A1 }$$

where η is the symmetric mass ratio, m_tot is the total mass of the two objects, and v_rel is the relative speed of the two objects at infinity. We only consider GW capture here since, for compact objects, the tidal capture cross section is much smaller than the GW capture cross section.

We can then calculate the cross section, using the standard formula for gravitational focusing

$$\begin{eqnarray} \sigma _{cs} &=& \pi r_{p, \max }^2 \left[1 + \frac{2Gm_{{\rm tot}}}{r_{p, \max } v_{\rm rel}^2} \right] \simeq \frac{2\pi G m_{{\rm tot}} r_{p,{\rm max}}}{v_{\rm rel}^2} \nonumber \\ &=& 1.3 \times 10^{23}\ {\rm cm}^2 (\eta /0.25)^{2/7} (m_{\rm tot}/20 \ M_\odot)^2\nonumber \\ && \times\,(v_{\rm rel}/84\, {\rm km\ s}^{-1})^{-18/7}, \end{eqnarray} \tag{ A2 }$$

where we have scaled this to fiducial values that will be used later. The rate of eccentric captures for a single galaxy with a central supermassive black hole (SMBH) of mass M_S and velocity dispersion σ_disp(M_S) is then

$$\begin{equation} \Gamma _{\rm gal}(M_S) \simeq \int _{r_{\rm min}}^{r_i} dr 4\pi r^2 n_1(r, M_S) n_2(r, M_S) \sigma _{cs} v_{\rm rel}, \end{equation} \tag{ A3 }$$

where n₁ and n₂ are the number densities of each type of object as a function of radius, $r_i \equiv {GM}_S{/}\sigma _{\rm disp}^2$ is the radius of influence of the BH, and r_min is approximately the radius inside which there is only a single object.

We take the density of objects to be the same as in an isothermal distribution, and the velocity distribution to be Maxwellian at each radius with relative speed $\langle v_{\rm rel}^2\rangle = 2 v_{\rm circ}^2 \equiv 2{GM}_S/r = 2\sigma _{\rm disp}^2 (r/r_i)^{-1}$, which captures the behavior well near the central BH, where the contribution to the rate is the largest. The number densities are then given by

$$\begin{equation} n_i(r) = \frac{N_i}{4 \pi r_{\rm dyn}^3}\left(\frac{r}{r_{\rm dyn}}\right)^{-2}, \end{equation} \tag{ A4 }$$

where r_dyn defines the dynamical radius inside which twice the mass of the central BH is contained such that N_i ≡ 2κ_iM_S/m_i and κ_i are the number of objects and total mass fraction of type i = 1, 2 within r_dyn, respectively. We take these scalings of n(r) and v_rel(r) for simplicity to highlight the sources of uncertainty and provide a rate estimate for the simplest case. For how different and more realistic scalings may alter the basic single galaxy rate, see the detailed discussion in O'Leary et al. (2009) and Appendix C of Kocsis & Levin (2012).

It is also convenient to define the geometric means of the number $\tilde{N}$ and mass fraction $\tilde{\kappa }$, such that

$$\begin{equation} \tilde{N}^2 \equiv N_1 N_2 = \frac{4 \tilde{\kappa }^2 M_S^2}{\eta \,m_{\rm tot}^2}, \quad \tilde{\kappa }^2 \equiv \kappa _1 \kappa _2, \end{equation} \tag{ A5 }$$

as well as the fiducial scalings σ₈₄ ≡ σ_disp/(84 km s⁻¹), η_0.25 ≡ η/0.25, and m₂₀ ≡ m_tot/(20 M_☉). Approximating $r_{\rm min} \simeq \tilde{N}^{-1} r_{\rm dyn}$, we find

$$\begin{eqnarray} &&\Gamma _{\rm gal}(M_S) \simeq \int _{r_{\rm min}}^{r_i} dr 4\pi r^2 \frac{\tilde{N}^2}{\left(4\pi r_{\rm dyn}^3\right)^2}\left(\frac{r}{r_{\rm dyn}}\right)^{-4} \sigma _{cs} v_{\rm rel}\qquad \end{eqnarray} \tag{ A6 }$$

$$\begin{eqnarray} \qquad &\simeq & 1.3 \times 10^{30}\ {\rm cm}^3 {\rm s}^{-1} \eta _{0.25}^{2/7} m_{20}^2 \sigma _{84}^{-11/7} \left(\frac{r_i}{r_{\rm dyn}}\right)^{-11/14}\nonumber\\ &&\times\, \frac{\tilde{N}^2}{4\pi r_{\rm dyn}^3} \int _{\tilde{N}^{-1}}^{r_i/r_{\rm dyn}} dx x^{-2 + 11/14}, \end{eqnarray} \tag{ A7 }$$

where we have averaged over the Maxwellian distribution, $\langle v_{\rm rel}(r_i)^{-11/7} \rangle = 1.15 \times \sigma _{\rm disp}^{-11/7}$, and changed integration variables to x ≡ r/r_dyn. Assuming that $(\tilde{N}r_i/r_{\rm dyn})^{3/14} \gg 1$, we have

$$\begin{eqnarray} \Gamma _{\rm gal}(M_S) &=& 1.2\times 10^{-10}\ {\rm yr}^{-1} (r_i/r_{\rm dyn})^{31/14}\nonumber\\ &&\times\,\eta _{0.25}^{-23/28} \tilde{\kappa }_{2.5}^{31/14} m_{20}^{-3/14} M_{4e6}^{-11/14} \sigma _{84}^{31/7}, \end{eqnarray} \tag{ A8 }$$

where M_4e6 ≡ M_S/(4 × 10⁶ M_☉) is scaled to the Milky Way, and $\tilde{\kappa }_{2.5} \equiv \tilde{\kappa }/(2.5\%)$ as in Kocsis & Levin (2012). Applying the M–σ relation, $M_{4e6} =\sigma _{84}^4$ (Tremaine et al. 2002), we finally have

$$\begin{eqnarray} \Gamma _{\rm gal}(M_S) &=& 1.2\times 10^{-10}\ {\rm yr}^{-1} (r_i/r_{\rm dyn})^{31/14}\nonumber\\ &&\times\,\eta _{0.25}^{-23/28} \tilde{\kappa }_{2.5}^{31/14} m_{20}^{-3/14} M_{4e6}^{9/28}. \end{eqnarray} \tag{ A9 }$$

Thus far, this agrees relatively well with Kocsis & Levin (2012) and O'Leary et al. (2009).

From the scatter in the inferred nuclear star cluster relaxation time T_r for the low σ_disp galaxies given in Figure 1 of Merritt et al. (2007), O'Leary et al. (2009) and Kocsis & Levin (2012) claim that the variance of the central number density scales the average rate per galaxy by a factor $\xi = \overline{n^2}/\overline{n}^2 \gtrsim 30$, increasing their total rate substantially.

Here, in calculating the average rate over many galaxies, we will carefully consider the effect of variation in both the M–σ relation and the scaling of the central density implied by Merritt et al. (2007)—related to the parameter r_i/r_dyn—and show that such a substantial increase in the inferred rate is not warranted.

Variation in the M–σ relation. We begin by taking the generous assumption that the intrinsic scatter in the M–σ relation is ∼0.5 dex in M_S (Tremaine et al. 2002). We take the distribution of M_S for a fixed σ_disp to be log-normal such that $M_{4e6} = C_{M\sigma }\times \sigma _{84}^4$, where C_Mσ is a random variable with log-normal probability distribution with geometric mean 〈C_Mσ〉 = 1 and scale factor $\delta _{M\sigma } = \ln \sqrt{10}$.

Our single galaxy rate (Equation (A8)) has scaling such that for a fixed SMBH mass bin, $\Gamma _{\rm gal}(M_S, C_{M\sigma }) \propto C_{M\sigma }^{-31/28}$. Using the properties of log-normal random variables, the scaling to the average rate over this distribution is given as

$$\begin{equation} \xi _{M\sigma } \equiv \frac{\overline{C_{M\sigma }^{-31/28}}}{\langle C_{M\sigma }\rangle ^{-31/28}} = \exp \left[ \frac{1}{2}\left(\frac{31}{28}\right)^2\delta _{M\sigma }^2\right] \simeq 2.25, \end{equation} \tag{ A10 }$$

where $\overline{f}$ denotes averaging of f over the distribution.

Variation in r_i/r_dyn. The ratio of the radius of influence, r_i, to the dynamical radius, r_dyn, determines the relative density of the nuclear cluster, and varies for different N-body models from ∼0.1 to 1. (Binney & Tremaine 2008). For simplicity we will take this ratio to be a log-normal distributed random variable independent of C_Mσ.

Figure 1 of Merritt et al. (2007) showed an estimate of the relaxation time at the dynamical radius as a function of the velocity dispersion σ_disp for nuclei of early-type galaxies in the ACS Virgo Cluster Survey (Côté et al. 2004). While the majority of the scatter in this distribution is for low-luminosity unresolved nuclear star clusters, an indication that much of this scatter may be due to observational uncertainty, for the purposes of this discussion we will assume, as in O'Leary et al. (2009), that this scatter is intrinsic. We will again assume, for simplicity, that the relaxation time T_r for fixed σ_disp is distributed log-normally, with a generous estimate of 1.5 dex standard deviation for the scatter in log T_r, such that the scale factor is $\delta _{T_r} \simeq 1.5 \ln 10$.

The nuclear relaxation time is given by

$$\begin{equation} T_r(r_{\rm dyn}) \simeq \frac{0.34 \sigma _{\rm disp}^3}{G^2 n(r_{\rm dyn}) \ln \Lambda } = \frac{0.34\sigma _{\rm disp}^3}{G^2 \ln \Lambda }\frac{4\pi r_{\rm dyn}^3}{2 \kappa M_S} \end{equation} \tag{ A11 }$$

(Spitzer 1987), where ln Λ is the Coulomb logarithm. For fixed σ_disp we can rewrite the T_r in terms of the random variables C_Mσ and r_i/r_dyn:

$$\begin{equation} T_r(r_{\rm dyn}) \simeq \frac{0.68\pi G M_S^2}{\ln \Lambda \kappa \sigma _{\rm disp}^3} \left(\frac{r_i}{r_{\rm dyn}}\right)^{-3} \sim \sigma _{\rm disp}^5 C_{M\sigma }^2 \left(\frac{r_i}{r_{\rm dyn}}\right)^{-3}. \end{equation} \tag{ A12 }$$

For independent log-normal random variables C_Mσ and r_i/r_dyn, this gives

$$\begin{equation} \delta _{T_r}^2 = (2\delta _{M\sigma })^2 + (3 \delta _{\rm dyn})^2, \end{equation} \tag{ A13 }$$

where δ_dyn is the standard deviation of ln r_i/r_dyn. This then gives the scaling due to variation in r_i/r_dyn to the average rate of

$$\begin{eqnarray} \xi _{\rm dyn}& \equiv & \frac{\overline{(r_i/r_{\rm dyn})^{31/14}}}{\langle r_i/r_{\rm dyn} \rangle ^{31/14}} = \exp \left[\frac{1}{2} \left(\frac{31}{14}\right)^2 \delta _{\rm dyn}^2 \right]\nonumber\\ & = & \exp \left[\frac{1}{2} \left(\frac{31}{14}\right)^2 \frac{\delta _{T_r}^2 - 4\delta _{M\sigma }^2}{9} \right] \simeq 6.1. \end{eqnarray} \tag{ A14 }$$

Final rate. The central BH mass function can be estimated by

$$\begin{eqnarray} \Phi (M_S) &\equiv& \frac{dn_{{\rm gal}}}{d \ln M_S} \simeq 0.0077\, {\rm Mpc}^{-3}\nonumber\\ && \times \left(\frac{M_S}{M_*}\right)^{\alpha + 1} \exp \left[ - \left(\frac{M_S}{M_*}\right)^{-\beta }\right] \end{eqnarray} \tag{ A15 }$$

(Shankar et al. 2004), where α ≃ −1.11, β ≃ 0.5, and M_* ≃ 6.4 × 10⁷ M_☉, assuming the local Hubble constant H_o = 70 km s⁻¹ Mpc⁻¹.

Integrating our mass dependence $M_{4e6}^{9/28}$ over the mass function, we can obtain the effective density

$$\begin{equation} n_{\rm gal, eff} = \int _{M_{\rm S, min}}^{M_{\rm S, max}} M_{4e6}^{9/28} \Phi (M_S) \frac{d M_S}{M_S} \end{equation} \tag{ A16 }$$

with which to multiply our single galaxy rate evaluated at M_4e6 = 1. The mass function (Equation (A15)) is only constructed to be valid between 10⁶ M_☉ ≲ M_S ≲ 5 × 10⁹ M_☉, however, we expect significant contribution from smaller galaxies. The halo mass function increases for lower mass, however, for dwarf galaxies the stellar mass (and therefore SMBH mass) to halo mass ratio drops significantly and the expected number density at that mass should also fall. If we integrate down to only M_{S, min} ∼ 10⁶ M_☉, then this gives us an effective density n_{gal, eff} ≃ 0.043 Mpc⁻³, while extending this mass function down to a cutoff of M_{S, min} ∼ 10⁴ M_☉ yields n_{gal, eff} ≃ 0.067 Mpc⁻³. With this in mind, we take the fiducial value of the effective density to be n_{gal, eff} = n_{gal, 5} × 0.05 Mpc⁻³.

Scaling to the Milky Way where O'Leary et al. (2009) assume a fiducial value of r_i/r_dyn ≃ 0.5, we have our final rate of eccentric captures in galactic nuclei,

$$\begin{eqnarray} &&\Gamma ^{({\rm GN}, {\rm EC})}_{\rm tot} = \frac{4}{3}\pi d^3 n_{\rm gal, eff} \overline{\Gamma _{{\rm gal}}}(M_{4e6} = 1)\qquad \end{eqnarray} \tag{ A17 }$$

$$\begin{eqnarray} \ \ \, \, \qquad &\simeq & 0.6 \,{\rm yr}^{-1}\left(\frac{\xi _{M\sigma }}{2.25}\right) \left(\frac{\xi _{\rm dyn}}{6.1}\right)\left(\frac{ \langle r_i/r_{\rm dyn}\rangle }{0.5}\right)^{31/14} \left(\frac{H}{r}\right)^{-2}\nonumber\\ && \times\,\eta _{0.25}^{-23/28} \tilde{\kappa }_{2.5}^{31/14} m_{20}^{-3/14} n_{\rm gal, 5} d_{2\,{\rm Gpc}}^3, \end{eqnarray} \tag{ A18 }$$

within d_2 Gpc × 2 Gpc, where we take 2 Gpc for the fiducial value as it is roughly the advanced LIGO horizon distance for 10 M_☉ BH–BH eccentric captures. Here, we have also included an additional factor (H/r)⁻², which may increase the density for nuclear clusters where significant flattening has occurred (B. Kocsis 2013, private communication). Assuming the generous fiducial values for ξ_Mσ and ξ_dyn, and no significant flattening, the rate for eccentric capture of 10 M_☉ BH–BH encounters is

$$\begin{equation} \Gamma _{\rm BH-BH}^{({\rm GN}, {\rm EC})} \simeq 0.02 \, {\rm yr}^{-1}\ {\rm Gpc}^{-3}. \end{equation} \tag{ A19 }$$

Less generous estimates for the intrinsic scatter, δ_Mσ ≃ 0.3ln 10 (Tremaine et al. 2002) and $\delta _{T_r} \simeq \ln 10$ give ξ_Mσ ≃ 1.33 and ξ_dyn = 2.52, which reduces the above rate by a factor of ∼4.

We are now (finally) ready to estimate the rates for NS eccentric captures in galactic nuclei. There is large uncertainty in the NS production rate in nuclear star clusters, mostly due to two factors. First, the IMF is unknown and could range from the standard Salpeter IMF to an extremely top-heavy IMF (e.g., Bartko et al. 2010). Second, there is great uncertainty in the effects of mass-loss for determining the fraction of stars with M_ZAMS ≳ 8 M_☉ that will become NSs, and the fraction that will collapse to form BHs (see, e.g., O'Connor & Ott 2011 for discussion).

With these considerations, we will assume that between ∼1% and 10% of the stars in a nuclear star cluster will become NSs. The low end of this range is for a standard IMF, with low mass loss and metallicity, while the high end roughly corresponds to a top-heavy IMF with more mass loss in the NS progenitors. Assuming an average stellar mass in nuclear star clusters of ∼0.5 M_☉, this gives κ_NS ∼ 0.03–0.3.

NSs have only had time to segregate in galaxies for which σ_disp ≲ 50 km s⁻¹ (Miller & Lauburg 2009). Unfortunately, these low mass galaxies also have escape velocities ∼2σ_disp, and thus NS kick velocities of ∼100 km s⁻¹ will significantly reduce the retention fraction after NS formation, much like in globular clusters. So we will assume the scaling above with no significant enhancements due to segregation or flattening.

Substituting m_tot = 2.8 M_☉, η = 0.25, and κ_NS = 0.03–0.3 into (A18) and continuing to use the very generous assumptions above for ξ_Mσ ≃ 2.25 and ξ_dyn ≃ 6.1, we get a rate for eccentric capture of

$$\begin{equation} \Gamma ^{({\rm GN}, {\rm EC})}_{\rm NS-NS} \simeq 0.04\hbox{--}6 \,{\rm yr}^{-1}\ {\rm Gpc}^{-3}. \end{equation} \tag{ A20 }$$

For 10 M_☉ BH–1.4 M_☉ NS encounters we take m_tot = 11.4 M_☉, η ≃ 0.11, κ_BH ≃ 0.025, and κ_NS ≃ 0.03–0.3, and get a rate

$$\begin{equation} \Gamma ^{({\rm GN}, {\rm EC})}_{\rm BH-NS} \simeq 0.05\hbox{--}0.6\, {\rm yr}^{-1}\ {\rm Gpc}^{-3}. \end{equation} \tag{ A21 }$$

However, in systems where the BHs dominate the core, a cusp of more massive objects tends to flatten out the distribution of lighter objects compared to an isothermal density profile. O'Leary et al. (2009) find that the distribution can approach a power law index of 1.5 for the lighter objects in their Fokker–Planck calculation. Repeating the above calculations for n_NS∝r^−3/2 as a lower bound, we find

$$\begin{equation} \Gamma ^{({\rm GN},{\rm EC})}_{\rm NS-NS} \simeq 0.003\hbox{--}0.3 \,{\rm yr}^{-1}\ {\rm Gpc}^{-3}, \end{equation} \tag{ A22 }$$

noting that the rate in a single galaxy is no longer dominated by the contribution due to the innermost objects.