Quick recipes for gravitational-wave selection effects

Let us first consider the case of a single detector. The optimal SNR is given by

$$\begin{equation} \rho_{\mathrm{opt}} = 2 \sqrt{\int_0^\infty \frac{|h\left(f\,\right)|^2}{S\left(f\,\right)} \mathrm{d} f} \,. \end{equation} \tag{ 1 }$$

where f indicates frequency, h(f) is the GW strain, and S(f) is the one-sided power spectral density of the detector. The word 'optimal' has sometimes been used in the literature (including by some of the authors [12]) to indicate the SNR of an optimally oriented source. In this paper, we refer to the more common meaning of optimality with respect to noise realizations.

A common approximation [8] is to take $\rho_{\mathrm{opt}}$ as a detection statistics, i.e. assume that a source is detectable if the optimal SNR is greater than a given threshold $\rho_{\mathrm{t}}$. That is, we write

$$\begin{equation} p_{\mathrm{det}}\left(\theta,\xi\right) = \mathcal{I}\left[\rho_{\mathrm{opt}} \left(\theta,\xi\right) > \rho_{\mathrm{t}}\right]\,, \end{equation} \tag{ 2 }$$

where $\mathcal{I}$ is an indicator function equal to one if the condition inside the brackets is true and zero otherwise ⁴ . For a given distribution of extrinsic parameters $p(\xi)$, one can compute (with an abuse of notation)

$$\begin{equation} p_{\mathrm{det}}\left(\theta\right) = \int p_{\mathrm{det}}\left(\theta, \xi\right) p\left(\xi\right) \mathrm{d}\xi \,. \end{equation} \tag{ 3 }$$

For a single detector, the optimal SNR factorizes as follows [8, 12]

$$\begin{equation} \rho_{\mathrm{opt}} \left(\theta,\xi\right) = \omega\left(\xi\right) \rho_{\mathrm{max}}\left(\theta\right)\,, \end{equation} \tag{ 4 }$$

where $0\unicode{x2A7D} \omega\unicode{x2A7D} 1$ is a projection parameter and $\rho_{\mathrm{max}}(\theta)$ is the optimal SNR of an optimally oriented source (i.e. a binary with face-on inclination located overhead the detector). In particular, one has $\xi = \{\iota,\vartheta,\phi,\psi\}$ and

$$\begin{align} \omega \, & = \; \Bigg\{\left(\frac{1+\cos^2\iota}{2} \right)^2 \left[ \frac{1}{2}\left(1+\cos^2\vartheta \right)\cos 2\phi \cos 2\psi - \cos \vartheta\sin 2\phi \sin 2\psi\right]^2 \nonumber \\ & \quad + \cos^2\iota \left[ \frac{1}{2}\left(1+\cos^2\vartheta \right) \cos 2\phi \sin 2\psi + \cos\vartheta\sin 2\phi \cos 2\psi \right]^2 \Bigg\}^{1/2} \,, \end{align} \tag{ 5 }$$

where ι is the binary inclination, ϑ and φ are a polar and an azimuthal angle for the sky location, and ψ is the polarization angle. Note that the factorization of equation (4) breaks down for precessing sources because the inclination of the orbital plane depends on the emitted frequency. That said, this effect is likely to be mild because current GW observations cover ${\lesssim} 1$ precession cycle.

With this factorization, the integral in equation (3) becomes [8]

$$\begin{equation} p_{\mathrm{det}}\left(\theta\right) = \!\int\! \mathcal{I}\left[\omega \rho_{\mathrm{max}} \left(\theta,\xi\right) > \rho_{\mathrm{t}}\right] p\left(\omega\right) \mathrm{d}\omega = \! \int_{\rho_{\mathrm{t}}/ \rho_{\mathrm{max}}\left(\theta\right)}^\infty\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!p\left(\omega\right) \mathrm{d}\omega\,. \end{equation} \tag{ 6 }$$

Assuming that the distribution $p(\xi)$ of the extrinsic parameters is known, this implies that selection effects can be estimated by evaluating the complementary cumulative distribution function of ω [8]

$$\begin{equation} p\left(\omega>\omega_0\right) = \int_{\omega_0}^\infty p\left(\omega^{^{\prime}}\right) \mathrm{d}\omega^{^{\prime}}\,. \end{equation} \tag{ 7 }$$

The simplicity of this approximation is appealing: estimating selection effects reduces to evaluating a single waveform h(f) for an optimally oriented source with intrinsic parameters θ, calculating $\rho_{\mathrm{max}}$ from equation (1), and evaluating the universal function $p({\gt}\omega)$ at $\omega = \rho_{\mathrm{t}}/ \rho_{\mathrm{max}}$.

Figure 1 shows the outcome of this procedure assuming that sources are distributed isotropically in the sky (i.e. we take uniform distributions in $\cos\iota$, $\cos\vartheta$, φ, and ψ). To facilitate comparison with the rest of the paper, we show the averaged detectability $p_{\mathrm{det}}(\theta)$ as a function of $\rho_{\mathrm{max}}/ \rho_{\mathrm{t}}$ instead of its inverse, even though the latter enters equation (6) directly. For instance, the black curve in the top panel of figure 1 indicates that at least ${\sim} 80\%$ of the binaries with a given set of intrinsic parameters θ will be detectable if at least one of these sources has an SNR that is ${\sim} 5$ times above the detection threshold.

Zoom In Zoom Out Reset image size

Because noise realizations, any specific GW source will not be observed with SNR $\rho_{\mathrm{opt}}$ but rather with some other value $\rho_{\mathrm{obs}}$. In the standard matched-filtering approach to GW detection, $\rho_{\mathrm{obs}}$ is given by the filtered signal (made of both GWs and noise) normalized by its own root-mean-square; see [21, 22]. One can thus include the effect of noise realizations in the estimate of $p_{\mathrm{det}}$ by thresholding the observed SNR $\rho_{\mathrm{obs}}$ instead of $\rho_{\mathrm{opt}}$ and computing the expectation value over noise realizations n, i.e.

$$\begin{equation} p_{\mathrm{det}}\left(\theta,\xi\right) = \int \mathcal{I}\left[\rho_{\mathrm{obs}}\left(n, \theta,\xi\right) > \rho_{\mathrm{t}}\right] \, p\left(n\right)\, \mathrm{d} n\,. \end{equation} \tag{ 8 }$$

Assuming the noise in the detector is Gaussian and stationary, the observed SNR is distributed normally around the optimal SNR with unit variance (see [21, 22] but also [10] for caveats). Because of this property, computing $\rho_{\mathrm{opt}}$ reduces to evaluating in equation (1) and adding a variance of one. One has

$$\begin{equation} p\left(\rho_{\mathrm{obs}}\right) = \frac{1}{\sqrt{2\pi}} \exp\left[ - \frac{\left(\rho_{\mathrm{obs}} - \rho_{\mathrm{opt}}\right)^2}{2} \right]\,, \end{equation} \tag{ 9 }$$

and thus [23]

$$\begin{align} p_{\mathrm{det}}\left(\theta,\xi\right) & = \frac{1}{\sqrt{2\pi}} \int_{\rho_{\mathrm{t}}}^{\infty} \exp\left\{- \frac{\left[\rho_{\mathrm{obs}} - \rho_{\mathrm{opt}}\left( \theta,\xi\right)\right]^2}{2} \right\} \mathrm{d}\rho_{\mathrm{obs}} \nonumber \\ & = \frac{1}{2}\left\{1 + \mathrm{erf}\left[\frac{\rho_{\mathrm{opt}}\left(\theta,\xi\right) - \rho_{\mathrm{t}}}{\sqrt{2}} \right] \right\}\,, \end{align} \tag{ 10 }$$

where $\mathrm{erf}$ is the error function. Marginalizing over the extrinsic parameters as above yields

$$\begin{align} p_{\mathrm{det}}\left(\theta\right) = \frac{1}{2}\left\{1 + \int \mathrm{erf}\left[\frac{\omega\rho_{\mathrm{max}}\left(\theta\right) - \rho_{\mathrm{t}}}{\sqrt{2}} \right] p\left(\omega\right) \mathrm{d}\omega \right\}\,. \end{align} \tag{ 11 }$$

Our results are shown in figure 1 assuming isotropic sources. Note that, unlike equation (6), thresholding the observed SNR does not result in a universal function of $\rho_{\mathrm{max}}(\theta)/\rho_{\mathrm{t}}$ but rather a one-parameter family of functions, where the additional parameter is ρ_t itself.

In practice, introducing the SNR variance due to noise realizations causes variations in $p_{\mathrm{det}}$ that are of $\mathcal{O}(1\%)$. The effect decreases as $\rho_{\mathrm{t}}$ increases: if the threshold is large, it is less likely that introducing a variance of one might turn a detectable event into a non-detectable event, or vice versa. Broadly speaking, we find that the GW detectability computed by thresholding the matched filter SNR is larger (smaller) than that obtained by thresholding the optimal SNR when signals are weak (loud). From figure 1, the transition between these two behaviors takes place at $\rho_{\mathrm{max}}\simeq 4 \rho_{\mathrm{t}}$. The largest deviations are found at SNRs $\rho_{\mathrm{max}}\simeq 2 \rho_{\mathrm{t}}$.

The luminosity distance d_L of astrophysical objects in the nearby Universe is distributed geometrically, $p(d_L)\propto d_L^{-2}$, which implies that $\rho_{\mathrm{max}}\propto 1/d_L$ is distributed as $p(\rho_{\mathrm{max}}) \propto {\rho_{\mathrm{max}}^{-4}}$ [24] (but note this will not be true for third-generation detectors [25]). This implies that the left part of the curve in figure 1 where $\rho_{\mathrm{t}}\simeq \rho_{\mathrm{max}}$ has a disproportionally larger impact on the overall population of detected sources. Therefore, thresholding the optimal SNR instead of the observed SNR has the net effect of underestimating $p_{\mathrm{det}}$ (i.e. the residuals in figure 1 are positive) and thus overestimating the astrophysical merger rate.

One caveat here is that we did not truncate the Gaussian distribution in equation (9) to impose $\rho_{\mathrm{obs}}\unicode{x2A7E} 0$. This is appropriate as long as the threshold value is much greater than the SNR variance, i.e. $\rho_{\mathrm{t}}\gg 1$.

Let us now consider a network of multiple detectors. The network SNR is the root sum square of the individual SNRs, i.e.

$$\begin{equation} \rho_{\mathrm{obs}} \left(\left\{n\right\}, \theta,\xi\right) = \sqrt{\sum_i^{N} \rho^2_{\mathrm{obs,i}} \left(n_i, \theta,\xi\right) }\,, \end{equation} \tag{ 12 }$$

where N is the number of interferometers. Each of the $\rho_{\mathrm{obs,i}}$'s is distributed normally around optimal values $\rho_{\mathrm{opt, i}}$ with unit variance. The square root of the sum of Gaussian variates with different means and unit variances is distributed according to the non-central χ distribution (which is also known as the generalized Rayleigh distribution) [26]; see also [10]. The probability density function of $\rho_{\mathrm{obs}}$ is

$$\begin{align} p\left(\rho_{\mathrm{obs}} \right) & = {\rho}_{\mathrm{opt}} \left( \frac{\rho_{\mathrm{obs}}}{{\rho}_{\mathrm{opt}}} \right)^{ N/2} \exp\left(-\frac{\rho_{\mathrm{obs}}^2 + {\rho}_{\mathrm{opt}}^2}{2}\right) \, I_{N/2 -1}\left({\rho}_{\mathrm{opt}} \rho_{\mathrm{obs}} \right) \,, \end{align} \tag{ 13 }$$

where $I_\nu(z)$ is the modified Bessel function of the first kind and order ν [27]. The parameter $\rho_{\mathrm{opt}}$ in equation (13) is root sum square of the individual optimal SNRs calculated as in equation (1), i.e.

$$\begin{equation} \rho_{\mathrm{opt}} \left(\theta,\xi\right) = \sqrt{\sum_i^{N} \rho^2_{\mathrm{opt,i}} \left(\theta,\xi\right) }\,. \end{equation} \tag{ 14 }$$

Note how the detection probability in equation (13) only depends on the combined quantity $\rho_{\mathrm{opt}}$ and not on how this SNR is distributed among the various instruments in the network.

We can now threshold the matched filter SNR $\rho_{\mathrm{obs}}$ and compute the expectation value over noise realizations as in equation (8). We find this can also be written down using special functions. In particular, one has

$$\begin{equation} p_{\mathrm{det}}\left(\theta,\xi\right) = Q_{N/2}\left[\rho_{\mathrm{opt}} \left(\theta,\xi\right), \rho_{\mathrm{t}} \right]\,, \end{equation} \tag{ 15 }$$

where

$$\begin{equation} Q_{\nu }\left(a,b\right) = a^{1-\nu}\int _{b}^{\infty }x^{\nu }\exp \left(-{\frac {x^{2}+a^{2}}{2}}\right)I_{\nu -1}\left(ax\right)\,\mathrm{d} x \end{equation} \tag{ 16 }$$

is the generalized Marcum Q-function of order ν [15–17], which is used in the field of digital communications (but see e.g. [28, 29] for some previous appearances in GW astronomy). In words, the generalized Marcum Q-function is the complementary cumulative distribution function of the non-central χ distribution. A two-line Python code to evaluate equation (16) is described in appendix and made available at [30].

In figure 2, we evaluate equation (15) for the case of N = 3 interferometers. As expected, thresholding the observed SNR smoothens the sharp step one would instead obtain when using the optimal SNR as detection statistics, allowing for some finite probability for events below (above) threshold to be detected (missed). It is interesting to note that this effect is not symmetric: binaries with optimal network SNR that is below threshold are more likely to be detected than binaries with optimal network SNR above threshold are to be missed. That is, thresholding the optimal SNR $\rho_{\mathrm{opt}}$ instead of the matched filter SNR $\rho_{\mathrm{obs}}$ underestimates $p_{\mathrm{det}}$, hence overestimates the intrinsic merger rate. This effect is enhanced by the expected astrophysical SNR probability $p(\rho_{\mathrm{opt}})\propto \rho_{\mathrm{opt}}^{-4}$ [24], which implies binaries are more likely to be found below than above threshold.

Zoom In Zoom Out Reset image size

We further quantify this detection/non-detection asymmetry as follows, borrowing terminology from that of a classification problem [31] where the actual outcome is given by $\rho_{\mathrm{obs}}\gt\rho_{\mathrm{t}}$ and the predicted outcome is given by the test $\rho_{\mathrm{opt}}\gt\rho_{\mathrm{t}}$. For a set of sources with SNRs distributed according to $p(\rho_{\mathrm{opt}})$, there are four mutually exclusive cases:

• $\rho_{\mathrm{t}}\lt\rho_{\mathrm{obs}}, \rho_{\mathrm{opt}}$ or 'true positives'. These events are flagged as detectable irrespectively of the thresholding strategy. The fraction of sources in this class is
  $$\begin{align} {\mathrm{TP}} & = \int Q_{N/2}\left(\rho_{\mathrm{opt}} , \rho_{\mathrm{t}}\right) \; \mathcal{I}\left(\rho_{\mathrm{opt}}>\rho_{\mathrm{t}}\right) \, p\left(\rho_{\mathrm{opt}}\right) \, \mathrm{d} \rho_{\mathrm{opt}} \end{align} \tag{ 17 }$$
• $\rho_{\mathrm{obs}}, \rho_{\mathrm{opt}}\lt\rho_{\mathrm{t}}$ or 'true negatives'. These events are flagged as non detectable irrespectively of the thresholding strategy. The fraction of sources in this class is
  $$\begin{align} {\mathrm{TN}} & = \int \left[1- Q_{N/2}\left(\rho_{\mathrm{opt}} , \rho_{\mathrm{t}}\right)\right]\; \mathcal{I}\left(\rho_{\mathrm{opt}} \lt \rho_{\mathrm{t}}\right) \, p\left(\rho_{\mathrm{opt}}\right) \, \mathrm{d} \rho_{\mathrm{opt}} \end{align} \tag{ 18 }$$
• $\rho_{\mathrm{opt}}\lt \rho_{\mathrm{t}}\lt\rho_{\mathrm{obs}}$ or 'false negatives'. These events are detectable but would be classified as non detectable if one neglects the SNR variance. The fraction of sources in this class is
  $$\begin{align} {\mathrm{FN}} & = \int Q_{N/2}\left(\rho_{\mathrm{opt}} , \rho_{\mathrm{t}}\right) \; \mathcal{I}\left(\rho_{\mathrm{opt}} \lt \rho_{\mathrm{t}}\right) \, p\left(\rho_{\mathrm{opt}}\right) \, \mathrm{d} \rho_{\mathrm{opt}} \end{align} \tag{ 19 }$$
• $\rho_{\mathrm{obs}}\lt \rho_{\mathrm{t}}\lt\rho_{\mathrm{opt}}$ or 'false positives'. These events are non detectable but would be classified as detectable if one neglects the SNR variance. The fraction of sources in this class is
  $$\begin{align} {\mathrm{FP}} & = \int \left[1- Q_{N/2}\left(\rho_{\mathrm{opt}} , \rho_{\mathrm{t}}\right) \right]\; \mathcal{I}\left(\rho_{\mathrm{opt}}>\rho_{\mathrm{t}}\right) \, p\left(\rho_{\mathrm{opt}}\right) \, \mathrm{d} \rho_{\mathrm{opt}} \end{align} \tag{ 20 }$$

Figure 3 shows the fractions FN and FP as a function of the threshold $\rho_{\mathrm{t}}$ and the number of detectors N. We consider a population with $\rho_{\mathrm{opt}}\in [1,100]$ distributed according to $p(\rho_{\mathrm{opt}})\propto \rho_{\mathrm{opt}}^{-4}$. Both fractions go to zero as $\rho_{\mathrm{t}}$ increase, corresponding to the curves of figure 2 approaching a step function: imposing a high detectability threshold makes it less likely for sources with a given optimal SNR to be scattered above or below threshold by a specific noise realization. The rate of false negatives is about an order of magnitude larger than that of false positives, confirming our point above. This difference increases with the number of detectors in the network, which is a consequence of equation (16).

Zoom In Zoom Out Reset image size

If needed, one can marginalize equation (15) over the extrinsic parameters as in equation (3). Note however that the factorization of equation (4) is not useful when considering multiple detectors because sources cannot be optimally oriented for all interferometers at the same time.

We estimate the impact of our findings on a toy population of GW events observable by the LIGO/Virgo network. We take a simple population $p_{\mathrm{pop}}(\theta, \xi | \lambda)$ where black-hole binaries have source-frame primary masses $m_1\in[5 M_\odot, 50 M_\odot]$ distributed according to $p(m_1)\propto m_1^{-2.3}$, source-frame secondary masses $m_2\in[5 M_\odot,m_1]$ distributed uniformly, redshifts $z\in[0,1]$ distributed uniformly in comoving volume and source-frame time $p(z)\propto (\mathrm{d} V_c / \mathrm{d} z)/ (1+z)$, spins magnitudes $\chi_{1,2}\in[0,1]$ distributed uniformly, and spin directions distributed isotropically. We assume a flat ΛCDM cosmological model with parameters from [32]. For the extrinsic parameters, we take uniform distributions in $\cos\iota$, $\cos\vartheta$, φ, and ψ as above. We consider a three-instrument network made of LIGO Livingston, LIGO Hanford, and Virgo at their design sensitivities [9] and compute optimal SNRs $\rho_{\mathrm{opt, i}}$ using the the IMRPhenomX approximant [33]. For this example, our threshold is set to $\rho_{\mathrm{t}} = 12$.

Figure 4 shows the resulting distributions of optimal network SNRs $\rho_{\mathrm{opt}}(\theta,\xi)$ and probabilities of detection $p_{\mathrm{det}}(\theta,\xi)$. A Monte Carlo estimate of the integral in equation (21) computed over the entire population returns $p_{\mathrm{det}}(\lambda)\sim 0.027$ (we used 10⁶ samples, so the error on this number is ∼10⁻³). If one only selects events with $|\rho_{\mathrm{opt}}(\theta,\xi) - \rho_{\mathrm{t}}|\lt1$ (i.e. those close to threshold), we find $p_{\mathrm{det}}(\lambda)\sim 0.496$. This should be compared with the analogous value ∼0.435 one would instead obtain with a simple cut on the optimal SNR. Marginalizing over noise realizations when estimating selection effects results in a higher $p_{\mathrm{det}}$ and mostly affects subpopulations of sources that are close to detection threshold.

Zoom In Zoom Out Reset image size

It is important to remember that the SNR (either optimal or observed) provides approximate information on the GW detectability. The quantity returned by GW detection pipelines is the false-alarm rate (FAR), which is indeed the statistics in used state-of-the-art population analyses [4] to both compile the list of events and estimate selection effects (other selection strategies are sometime used, see e.g. $p_{\mathrm{astro}}$ in [34]). Using the population average of equation (21), we now present a calibration of the SNR threshold that enters our $p_{\mathrm{det}}$ approximation to reproduce the response of current detection pipelines as a function of their FARs.

We use software injections performed in real noise from the third LIGO/Virgo/KAGRA observing run [35]. They report FAR values for five different detection pipelines and LIGO-only (i.e. N = 2) optimal SNRs for a fiducial population of sources (some of the injections also include Virgo; we neglect this difference and have verified it has a negligible impact on our results). Note their adopted population is different from that of our toy case above (see [35] for details). We use the injection dataset labeled as 'mixture,' which is their recommended default. We consider all the injections provided, including cases where the FAR could not be quantified (and is reported as $\infty$ in the dataset). Note the injections were performed uniformly in time, thus taking into account the duty cycle of the detectors.

We compute the population-average detection probability by thresholding their FARs and match it with the population-average detection probability obtained with our $p_{\mathrm{det}}$ approach. Figure 5 shows the resulting mapping between the two quantities. We repeat our calculation for each of the five detection pipelines provided in the available dataset [35] as well as by selecting the minimum FAR for each injected source. The latter is in line with the criterion used in [4] for selecting compact binaries of astrophysical origin.

Zoom In Zoom Out Reset image size

This calculation is not reliable for SNRs that are ${\lesssim} 6$ because those injections were deemed 'hopeless' to save computational time [35]. For those values of $\rho_{\mathrm{t}}$, there is a (potentially large) number of missing sources with optimal SNRs below threshold that could have been scattered above threshold by the SNR variance. We thus restrict our analysis to $\rho_{\mathrm{t}}\unicode{x2A7E} 8$, which is a few standard deviations away from the 'hopeless' cut. Additionally, some pipelines have minimum and maximum FAR values they attempt to quantify. We estimated these limits from the provided datasets and truncated the curves in figure 5 accordingly. These limitations do not significantly impact the minimum FAR calculation across pipelines, at least in the region of interest shown in figure 5.

Figure 5 shows that selection effects computed using a minimum FAR threshold of ${\sim} 1$ yr⁻¹ are reproduced by an SNR-based cut with threshold $\rho_{\mathrm{t}}\simeq 11$. We find once more that the bias induced by thresholding $\rho_{\mathrm{opt}}$ instead of $\rho_{\mathrm{obs}}$ is of a few % and underestimates $p_{\mathrm{det}}$ across the entire range: that is, the $\rho_{\mathrm{t}}$ threshold for a given FAR obtained when marginalizing over noise realization is larger compared to the case where noise is neglected. The difference between the two treatments is smaller than the typical difference between the various detection pipelines (and indeed current LIGO/Virgo selection criteria combine information from multiple pipelines). The population used in [35] is deliberately broad, while our results above indicate that the subpopulation of sources that are close to threshold will be affected more significantly by the SNR variance. Confirming this expectation with dedicated pipeline injections is left to future work.

As shown in figure 5, we empirically find that the mapping between $\rho_{\mathrm{t}}$ and the minimum FAR threshold across the available pipelines is well described by a power law:

$$\begin{align} \log_{10} \left(\frac{{\mathrm{FAR}}_{\mathrm{t}}}{{\mathrm{yr}}^{-1}}\right) = p_1\, \rho_{\mathrm{t}} +p_0. \end{align} \tag{ 22 }$$

A simple least-square regression returns $p_1 = -1.80$ and $p_0 = 20.1$ when thresholding using $\rho_{\mathrm{obs}}$ (i.e. this is a fit to the solid black curve in figure 5), and $p_1 = -1.72$ and $p_0 = 18.7$ when thresholding using $\rho_{\mathrm{opt}}$ (i.e. this is a fit to the dashed black curve in figure 5). These fits can be used to quickly filter synthetic distributions according to the desired purity of the resulting GW simulated catalog, though with the important caveat that this relationship was calibrated on a specific population of sources. We stress that figure 5 and equation (22) present a calibration between thresholds, not a mapping of SNR to FAR for a given trigger and pipeline.