Rapid Premerger Localization of Binary Neutron Stars in Third-generation Gravitational-wave Detectors

Since the first direct gravitational-wave (GW) detection of the coalescence of binary neutron star (BNS) GW170817(Abbott et al. 2017a) and its electromagnetic (EM) counterparts (Abbott et al. 2017b), multimessenger observation of coalescing compact binaries has become an important tool for astrophysics. Joint GW-EM observations can provide a comprehensive understanding of the formation and evolution of BNS, and shed light on physics around compact objects (Cowperthwaite et al. 2017; Kasen et al. 2017; Margalit & Metzger 2017; Nicholl et al. 2017; Soares-Santos et al. 2017; Abbott et al. 2018; Annala et al. 2018; De et al. 2018; Mooley et al. 2018; Abbott et al. 2019; Capano et al. 2020).

Rapid GW detection and accurate localization is key to joint GW-EM observations, as most EM facilities need the direction from the GW observation. In addition to capturing the afterglow of BNS coalescences, early EM observations could offer unique insights of phenomena in BNSs that happen prior to or near the merger (e.g., tidal disruptions, magnetosphere interactions, r-process nucleosynthesis) and help build a picture of the entire evolution of kilonova in multiple frequency bands (Metzger & Piro 2014; Nicholl et al. 2017; Radice et al. 2018; Metzger 2020; Most & Philippov 2020; Cooper et al. 2023). Early detection and localization of BNS are therefore of great importance in GW astronomy. It has previously been demonstrated that there is a nonzero probability of detecting and localizing premerger BNS events with current (Magee et al. 2021; Kovalam et al. 2022) and near-future (Nitz et al. 2020; Sachdev et al. 2020; Banerjee et al. 2023; Magee & Borhanian 2022) GW observatories, typically within seconds to 1 minute before merger. Chaudhary et al. (2023) has recently investigated early warning for the fourth observing run of LIGO-Virgo-KAGRA collaboration, with multiple detection pipelines already equipped for early-warning searches (Messick et al. 2017; Nitz et al. 2018; Aubin et al. 2021; Chu et al. 2022). Machine-learning-based methods for early-warning detection are also making rapid progresses (Baltus et al. 2021; Wei & Huerta2021; Yu et al. 2021).

Being limited by low sensitivities below 20 Hz, the BNS signal is only detectable for ∼1 minute in current GW detectors. Given the communication time delay between multimessenger community and the ∼10–100 s slew time of modern telescopes (Banerjee et al. 2023), it is basically impossible to capture premerger or near-merger transients from BNS coalescence without forewarning. Several third-generation (3G) GW detectors have been proposed, including Einstein Telescope (ET; Punturo et al. 2010) and Cosmic Explorer (CE; Reitze et al. 2019; Evans et al. 2023), with low-frequency sensitivities significantly improved. These would allow us to detect BNS signals more than 30 minutes before the merger, rendering precise early-warning localization possible (Chan et al. 2018; Akcay 2019; Nitz & Dal Canton 2021; Borhanian & Sathyaprakash 2022; Ronchini et al. 2022; Branchesi et al. 2023).

However, data analysis of BNS in 3G detectors can be challenging. The long signal makes matched filtering extremely expensive to perform, and it is modulated by changes of the antenna response functions due to the Earth's rotation. Neglecting the Earth's rotation should have little impact on detection, as signal-to-noise ratio (S/N) is mostly contributed by the last stages of the signal (see Figure 1). However, ignoring it causes loss of information and could lead to biases in parameter estimation, especially for sky location parameters. In this work, we demonstrate that multiband analysis (Cannon et al. 2012; Adams et al. 2016; Aubin et al. 2021) is an effective way of solving these issues, and fast localization algorithms can be built upon multiband detection statistics.

Zoom In Zoom Out Reset image size

Multiband analysis is based on the fact that orbital evolution of the quasi-circular BNS inspiral stage is well modeled. Observable BNSs are not likely to have large spins (Abbott et al. 2017a, 2019, 2020) or precession, and the frequency evolves as

$$\begin{eqnarray}&&f(\tau )=134\,\mathrm{Hz}{\left(\displaystyle \frac{1.21{M}_{\odot }}{{{ \mathcal M }}_{c}}\right)}^{5/8}{\left(\displaystyle \frac{1\,{\rm{s}}}{\tau }\right)}^{3/8}\end{eqnarray} \tag{ 1 }$$

to the Newtonian order (Maggiore 2007), where f is GW frequency, τ is the time before merger, and ${{ \mathcal M }}_{c}$ is the chirp mass of the binary. The monotonic evolution of GW frequency ensures a one-to-one correspondence between time segments and frequency bands, i.e., chopping the GW waveform into multiple time segments, [t_n , t_n−1), [t_n−1, t_n−2),...,[t₁, t₀), results in a corresponding sequence of frequency bands [f_n , f_n−1), [f_n−1, f_n−2),...,[f₁, f₀) defined via Equation (1). One can choose the length of time intervals such that within each interval the Earth's rotation can be ignored, i.e., the detectors' antenna response functions can be assumed constant and current matched filtering techniques can be directly employed. One can also down-sample the data in each frequency band according to its highest frequency to reduce computational cost. Results from each time segment can be combined in succession as new data come in.

Figure 1 shows the multiband scheme in this work. We consider the negative latency up to 60 minutes, and choose 2 minute segments and 256 Hz sampling frequency until the final 2 minutes. In the last 2 minutes the S/N grows rapidly, and the GW reaches high frequencies; therefore a finer time resolution is used to improve the detection and localization. In addition to the limit from Nyquist–Shannon sampling theorem, in practice we find that a sampling rate that is higher than Nyquist frequency could be helpful in localizing high-S/N events. As a demonstration, we equally divide the last 2 minutes, and employ sampling frequencies of 1024 and 4096 Hz, respectively. A more elaborate segmentation is also sensible, e.g., (Morisaki 2021), but we leave a comprehensive investigation of the multiband scheme to future works.

We will build a fast localization algorithm based on the above multibanding scheme. Given the large number of BNS detections in the 3G detectors (typically ∼10⁵ per year; Borhanian & Sathyaprakash 2022; Branchesi et al. 2023), an efficient and light algorithm is necessary. While machine-learning methods (Dax et al. 2021; Chatterjee & Wen 2022; Chatterjee et al. 2022; Gabbard et al. 2022) are gradually making progresses, Bayestar (Singer & Price2016), which performs a five-fold numerical marginalization over nuisance extrinsic parameters, has been used as the standard low-latency localization approach throughout the observing runs of the current generation of GW detectors. It takes Bayestar ∼2 s to generate a sky map with 32 threads (Singer & Price 2016), and lower latency can be achieved with more CPUs or narrower bandwidth. For instance, Magee et al. (2021) demonstrated ∼0.5 s of computation time of Bayestar with sufficient (>100) threads for early-warning triggers, and ∼1.1 s for full-bandwidth triggers. In this work, to reduce the computational cost of dealing with vast number of events, we will make use of a semianalytical localization algorithm for GWs (SealGW; Hu et al. 2021) that has recently been implemented in the SPIIR detection pipeline (Chu et al. 2022) and is publicly available. ¹ SealGW performs a semianalytical marginalization over nuisance parameters, achieves a faster performance than Bayestar, and retains a reasonable accuracy. A more detailed description will be given in the following sections.

As intrinsic parameters of GWs are initially estimated by matched filtering searches and their errors are semiindependent with errors in sky localization (Singer & Price 2016), they are treated as perfectly known in online fast localization. The Bayesian posterior probability distribution for extrinsic parameters can be written as

$$\begin{eqnarray}&&p({{\boldsymbol{\vartheta }}}_{\mathrm{ex}}| d,{{\boldsymbol{\vartheta }}}_{\mathrm{in}})\propto p({{\boldsymbol{\vartheta }}}_{\mathrm{ex}})p(d| {{\boldsymbol{\vartheta }}}_{\mathrm{ex}},{{\boldsymbol{\vartheta }}}_{\mathrm{in}}),\end{eqnarray} \tag{ 2 }$$

where d is the data, p( ϑ _ex) is the prior distribution, and p(d∣ ϑ _ex, ϑ _in) is the likelihood. Among ϑ _ex = {α, δ, t_c , ψ, r, ϕ_c , ι}, R.A. α and decl. δ describe the source sky direction while the other nuisance parameters should be marginalized. Direct analytical marginalization is impossible, but after a parameter conversion

$$\begin{eqnarray}&&{{\boldsymbol{\vartheta }}}_{\mathrm{ex}}=\{\alpha ,\delta ,{t}_{c},\psi ,r,{\phi }_{c},\iota \}\to {{\boldsymbol{\vartheta }}}_{\mathrm{ex}}^{{\prime} }=\{\alpha ,\delta ,{t}_{c},{\boldsymbol{A}}\},\end{eqnarray} \tag{ 3 }$$

where

$$\begin{eqnarray}\begin{array}{l}A=\left(\begin{array}{cc}\cos 2\psi & \sin 2\psi \\ -\sin 2\psi & \cos 2\psi \end{array}\right)\left(\begin{array}{cc}\displaystyle \frac{1+{\cos }^{2}\iota }{2r} & \\ & \displaystyle \frac{\cos \iota }{r}\end{array}\right)\\ \ \ \times \ \left(\begin{array}{cc}\cos {\phi }_{c} & \sin {\phi }_{c}\\ -\sin {\phi }_{c} & \cos {\phi }_{c}.\end{array}\right),\end{array}\end{eqnarray} \tag{ 4 }$$

the likelihood function appears as a Gaussian in the matrix A and is therefore analytically tractable (Hu et al. 2021). The posterior of the source sky location takes the form

$$\begin{eqnarray}&&\begin{array}{l}p(\alpha ,\delta | d,{{\boldsymbol{\vartheta }}}_{\mathrm{in}})=\int {{dt}}_{c}{d}^{4}{\boldsymbol{A}}p({{\boldsymbol{\vartheta }}}_{\mathrm{ex}}^{{\prime} }| d,{{\boldsymbol{\vartheta }}}_{\mathrm{in}})\\ \quad =\ \int {{dt}}_{c}I(\rho ({t}_{c}),\alpha ,\delta ),\end{array}\end{eqnarray} \tag{ 5 }$$

where ρ(t_c ) is the S/N time series. The analytical expression of I(ρ(t_c ), α, δ), including a comprehensive introduction and tests of SealGW, can be found in Hu et al. (2021).

Filtering each time series in the multiband scheme produces a separate complex S/N time series that must be coherently combined to achieve a precise localization (Allen et al. 2012; Singer & Price 2016). The gain in precision for long signals comes not only from the accumulation of absolute S/N, but also from phase drifts of the S/N time series due to the Earth's rotation because the rotation induces a longer equivalent network baseline (Wen & Chen 2010; Zhao & Wen 2018; Baral et al. 2023). A direct combination (linear addition) of multiple S/N time series is feasible for short signals, as it only requires a set of combination parameters to align the template bands in time and phase (Adams et al. 2016). However, for long signals, the direct combination scheme is no longer coherent because of the phase drifts due to the Earth's rotation. Including the phase drifts in combination parameters can lead to a coherent addition, but it loses information contained in the changing time delays between detectors, which depends on the sky location. Therefore, instead of directly adding S/N, we multiply likelihoods from every band before marginalization over nuisance extrinsic parameters. Since there is little overlap between bands (contributed by noise correlation and template overlaps, which are windowed out), they can be treated as independent measurements, as if there are different detectors at different frequency bands.

To assess the performance of the above localization scheme, we simulate a mock BNS catalog, assuming a three-detector network with one triangle ET at Virgo site and two L-shaped CEs at the LIGO Hanford and Livingston site, respectively. We use an analytical astrophysical population (Oguri 2018)

$$\begin{eqnarray}&&{R}_{\mathrm{obs}}(z)=\displaystyle \frac{{a}_{1}{e}^{{a}_{2}z}}{{e}^{{a}_{3}z}+{a}_{4}}\displaystyle \frac{1}{1+z}\displaystyle \frac{{{dV}}_{{\rm{c}}}}{{dz}}{\mathrm{Gpc}}^{-3}\,{\mathrm{yr}}^{-1}.\end{eqnarray} \tag{ 6 }$$

Here V_c is the comoving volume and we employ Planck15 cosmology (Ade et al. 2016). a_{2,3,4} are set to be {1.6, 2.1, 30} to model a peak at z ∼ 2. a₁ is scaled to match the local BNS merger rate given by Abbott et al. (2021), (${{ \mathcal R }}_{\mathrm{obs}}(z=0)\,={320}_{-240}^{+490}{\mathrm{Gpc}}^{-3}\,{\mathrm{yr}}^{-1}$). We use the estimated median value of the merger rate in this work. We simulate 68,000 BNS sources within z = 3, which corresponds to roughly one month of observations. We assume neutron star mass is uniformly distributed in [1.1M_⊙, 2M_⊙] in the source frame and isotropic sky distribution and inclination. Note that the location and configuration of detector networks are not settled yet and are subject to change, and the BNS mass distribution and the current merger rate density estimate have large uncertainties due to the as yet small number of BNS detections.

Signals are injected into Gaussian noise realizations and analyzed individually, i.e., we do not consider them to be overlapping with each other. Overlapping signals could cause dominant biases in parameter estimation, but mainly in the case when the merger times are very close (Himemoto et al. 2021; Relton & Raymond 2021; Samajdar et al. 2021; Pizzati et al. 2022; Relton et al. 2022; Hu & Veitch 2023a). Among the 68,000 BNS events evenly distributed in one month, roughly 1.3% of them have another event ending <0.5 s afterwards. Many of the signals are not actually detectable, further reducing the chance of significant bias. Even though the number of overlapping signals could be large given the large number of events expected to be detected with 3G detectors, our simulation will still apply to the vast majority of noninterfering signals.

We use the waveform model TaylorF2 (Blanchet et al. 1995; Poisson 1998; Buonanno et al. 2009) to generate GW signals and map frequency to time before merger via the stationary phase approximation with Equation (1). The 3.5 post-Newtonian waveform is a reasonable choice for analyzing quasi-circular inspiralling compact binaries (Faye et al. 2012). Several works suggest some non-quasi-circular binaries, like precessing and eccentric systems, or systems with strong higher-order emission, can be better localized (Ma et al. 2017; Kapadia et al. 2020; Singh et al. 2021; Tsutsui et al. 2021; McIsaac et al. 2023). However, that would require novel search algorithms (e.g., Fairhurst et al. 2020) upon which new fast localization methods would have to be built, because current localization methods, including Bayestar and SealGW, are based on aligned-spinning waveform templates in which plus and cross polarizations of GWs only have a phase difference.

We perform matched filtering assuming a perfect knowledge of intrinsic parameters and set total S/N > 12 as the detection criterion, where total S/N is converted from the multiband matched filtering outputs by the analytical expression of S/N at Newtonian order (Cannon et al. 2012). Matched filtering with known injection parameters is the ideal case, while in a realistic scenario one should build a template bank that achieves a reasonable match (e.g., >97%) everywhere in the parameter space. The purpose of this work is to assess the performance of the multiband localization scheme. We leave a dedicated long signal early-warning pipeline and simulations with more realistic mock data to future work.

For each simulation, we perform multiband matched filtering from 60 minutes before merger with low-frequency cutoff at 5 Hz. Figure 2 shows the cumulative number of events for different negative latencies. ∼10 events can be localized within 100 deg² 20 minutes before merger, and 6 minutes before merger the number of events increases to ∼100. Also, ∼1–10 events can be localized within 10 deg² up to 6 minutes before merger.

Zoom In Zoom Out Reset image size

Extreme early warnings are possible. Several events in our simulation are detected and preliminarily localized 40–50 minutes before merger and this number could be underestimated since our analysis has hard cutoffs at 5 Hz and 1 hr negative latency. ET would be able to collect sensible data down to ∼3 Hz and trigger even earlier detections (Borhanian & Sathyaprakash 2022; Branchesi et al. 2023). However, BNS with high negative latencies are not likely to be well localized until more data comes in, bringing higher S/Ns, wider frequency bands, and a longer equivalent network baseline. Multiband analysis helps update detection statistics and sky maps on the fly. Figure 3 shows the evolution of sky maps and localization areas in our simulation. The example sky map is from a 1.4 + 1.4M_⊙ BNS at 1000 Mpc detected 30 minutes before merger. It presents nested contours with new bands combined in succession and is finally pinpointed within 0.2 deg², but is already well localized ∼10 minutes before the merger. The localization area traces in the lower panel show the decreasing rate of localization areas: those localized within 100 deg² ∼20 minutes before merger in Figure 2 are generally detected 40–50 minutes before merger.

Zoom In Zoom Out Reset image size

Figure 4 is the P–P plot of our localization simulation, showing x% confidence region (x-axis) is able to include y% of total events (y-axis, scaled). The diagonal shapes shows the multiband localization scheme is reasonably self-consistent. The lines for 30+ minutes before merger have larger statistical fluctuations due to the insufficient number of samples.

Zoom In Zoom Out Reset image size

We tested the time cost of SealGW and Bayestar calculation with the same data (full-bandwidth ET+2CE network) and sky map resolution (n_side = 2048, finest pixel = 0.0008 deg²), as showed in Figure 5. Tests are performed on a 2.44 GHz processor with OpenMP multithreading. Thanks to the semianalytical property, SealGW can achieve ∼26 times faster speed than Bayestar with fewer threads, and the speed up factor goes down to ∼4 when more threads come in as the nonparallelizable calculation begins to dominate SealGW run time. It only takes SealGW ∼3 s with 1 thread and ∼0.5 s with 8 threads, which means SealGW is able to perform real-time localization with a low hardware requirement. Note that the time cost can be further reduced with narrower bandwidth or coarser sky map resolutions, e.g., time cost of SealGW can be halved when the finest resolution is 0.013 deg². The efficiency and cheapness is suitable for the 3G detector scenario in which number of detection can be huge. Nevertheless, a thorough estimate of early-warning latency would require a comprehensive design of detection pipeline structure, and there would be a wall time of ∼0.1 s to read and preprocess the data from pipeline outputs.

Zoom In Zoom Out Reset image size

We provide an exploratory demonstration of early-warning localization of long signals for 3G GW detector networks. We simulate a mock catalog for one month of observation with an ET+2CE network, and perform multiband analysis with the fast localization algorithm SealGW. We show that this is a efficient scheme for premerger localization.

Multiband analysis allows us to detect BNS in an early stage and update the results regularly with incoming data. There are tens of BNS detected more than 30 minutes before merger in our simulation, and localized within 100 deg² at ∼10 minutes before merger. 10 deg² can be achieved ∼6 minutes before merger. Since wide-field optical transient facilities usually have field of view of 1–10 deg² (see summaries in Sachdev et al. 2020; Ronchini et al. 2022), the precise premerger localization of BNS would be extremely helpful to finding EM counterparts before the merger and observing the entire process of BNS coalescence.

Our work here presents a solution for the crucial step of performing real-time localization in the context of online searches in 3G detectors, that effectively reduces the latency and computational burden arising from premerger localization. However, there remains the larger issue of developing the surrounding infrastructure to search for premerger signals and disseminating sky maps in low latency to observatories before detection. As an exploratory demonstration, we have made several simplifications to the problem, such as ignoring overlapping signals, assuming perfect matched filtering, and a relatively naive waveform segmentation. The merger rate estimation of BNS is also uncertain to date; therefore the absolute detection numbers should be interpreted as an order-of-magnitude estimation. We plan to explore the multiband analysis on a real detection pipeline and use a more accurate astrophysical population (which should be available in years with new observations) in our future work.

We thank Linqing Wen, Daniel Tang, and Chayan Chatterjee for helpful discussions. We are grateful for computational resources provided by Cardiff University, and funded by STFC grant ST/I006285/1. Q.H. is supported by CSC. J.V. is supported by STFC grant ST/V005634/1.

The sky map fits files for 30, 20, 10, 6, and 1 minute early-warning and postmerger triggers are openly available in zenodo at doi:10.5281/zenodo.8297806 (Hu & Veitch 2023b). Further data and scripts for reproducibility are available upon reasonable request.