Rapid pre-merger localization of binary neutron stars in third generation gravitational wave detectors

Pre-merger localization of binary neutron stars (BNSs) is one of the most important scientific goals for the third generation (3G) gravitational wave (GW) detectors. It will enable the electromagnetic observation of the whole process of BNS coalescence, especially for the pre-merger and merger phases which have not been observed yet, opening a window for deeper understandings of compact objects. To reach this goal, we describe a novel combination of multi-band matched filtering and semi-analytical localization algorithms to achieve early-warning localization of long BNS signals in 3G detectors. Using our method we are able to efficiently simulate one month of observations with a three-detector 3G network, and show that it is possible to provide accurate sky localizations more than 30 minutes before the merger. Our simulation shows that there could be ~ 10 (~ 100) BNS events localized within 100 deg2, 20 (6) minutes before merger, per month of observation.


INTRODUCTION
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 lights on physics around compact objects (Abbott & et al 2019;De et al. 2018;Mooley et al. 2018;Capano et al. 2020;Nicholl et al. 2017;Abbott et al. 2018;Annala et al. 2018;Kasen et al. 2017;Margalit & Metzger 2017;Cowperthwaite et al. 2017;Soares-Santos et al. 2017).
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 BNS 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 (Cooper et al. 2022;Most & Philippov q.hu.2@research.gla.ac.ukJohn.Veitch@glasgow.ac.uk 2020; Nicholl et al. 2017;Metzger & Piro 2014;Radice et al. 2018;Metzger 2020).Early detection and localization of BNS are therefore of great importance in GW astronomy.It has previously been demonstrated that there is a non-zero probability of detecting and localizing pre-merger BNS events with current (Kovalam et al. 2022;Magee et al. 2021) and near-future (Sachdev et al. 2020;Nitz et al. 2020;Magee & Borhanian 2022;Banerjee et al. 2022) GW observatories, typically within seconds to one 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 (Chu et al. 2022;Nitz et al. 2018;Messick et al. 2017;Aubin et al. 2021).Machine learning based methods for early warning detection are also making rapid progresses (Baltus et al. 2021b;Wei & Huerta 2021;Yu et al. 2021;Baltus et al. 2021a) 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 multi-messenger community and the ∼ 10 − 100 seconds slew time of modern telescopes (Banerjee et al. 2022), it is basically impossible to capture pre-merger or near-merger transients from BNS coalescence without fore-warning.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.
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 (SNR) is mostly contributed by the last stages of the signal (see Fig. 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 multi-band analysis (Aubin et al. 2021;Adams et al. 2016;Cannon et al. 2012) is an effective way of solving these issues, and fast localization algorithms can be built upon multi-band detection statistics.
Multi-band 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. 2020;Abbott & et al 2019Abbott & et al , 2017a) ) or precession, and the frequency evolves as (1) to the Newtonian order (Maggiore 2007), where f is GW frequency, τ is the time before merger, and 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 1 , t 0 ), results in a corresponding sequence of frequency bands 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 comes in.Fig. 1 shows the multi-band scheme in this work.We consider the negative latency up to 60 min, and choose 2 min segments and 256 Hz sampling frequency until the final two minutes.In the last two minutes the SNR grows rapidly, and the GW reaches high frequencies, therefore a finer time resolution is used to improve the An illustration of a chopped GW waveform that is alternately colored for different waveform bins and sampling frequencies.We use two-minute segments with 256 Hz sampling frequency for waveforms from 60 minutes before merger to two minutes before merger, and one-minute segments for the last two minutes with 1024 Hz and 4096 Hz sampling frequencies, respectively.Right axis: The cumulative SNR of the signal, showing the contribution of SNR from each time segment.
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 SNR events.As a demonstration, we equally divide the last two minutes, and employ sampling frequencies of 1024 Hz and 4096 Hz, respectively.A more elaborate segmentation is also sensible e.g.(Morisaki 2021), but we leave a comprehensive investigation of the multi-band scheme to future works.
We will build a fast localization algorithm based on the above multi-banding scheme.Given the large number of BNS detections in the 3G detectors (typically ∼ 10 5 per year (Branchesi et al. 2023;Borhanian & Sathyaprakash 2022)), an efficient and light algorithm is necessary.While machine learning methods (Gabbard et al. 2022;Dax et al. 2021;Chatterjee & Wen 2022;Chatterjee et al. 2022) are gradually making progresses, Bayestar (Singer & Price 2016), which performs a five-fold numerical marginalization over nuisance extrinsic parameters, has been used as the standard lowlatency localization approach throughout the observing runs of the current generation of GW detectors.It takes Bayestar ∼ 2 s to generate a skymap 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 semi-analytical 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 1 .SealGW performs a semi-analytical 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.

LOCALIZATION FOR LONG SIGNALS
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 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 , ι}, right ascension α and declination δ describes the source sky direction while the other nuisance parameters should be marginalized.Direct analytical marginalization is impossible, but after a parameter conversion where ) 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 where ρ(t c ) is the SNR timeseries.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 timeseries in the multi-band scheme produces a separate complex SNR 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 SNR, but also from phase drifts of the SNR timeseries due to the Earth's rotation because the rotation induces a longer equivalent network 1 https://git.ligo.org/spiir-group/SealGWbaseline (Wen & Chen 2010;Zhao & Wen 2018;Baral et al. 2023).A direct combination (linear addition) of multiple SNR timeseries 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 SNR, 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.

CATALOG SIMULATION
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 LIGO Hanford and Livingston site, respectively.We use an analytical astrophysical population (Oguri 2018) a 1 e a2z e a3z + a 4 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 1 is scaled to match the local BNS merger rate given by Abbott et al.
(2021) (R obs (z = 0) = 320 +490 −240 Gpc −3 yr −1 ).We use the estimated median value of the merger rate in this work.We simulate 68000 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 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;Pizzati et al. 2022;Relton & Raymond 2021;Relton et al. 2022;Samajdar et al. 2021;Hu & Veitch 2023a).Among the 68000 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 non-interfering 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 Eq. 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-quasicircular binaries, like precessing and eccentric systems, or systems with strong higher order emission, can be better localized (McIsaac et al. 2023;Tsutsui et al. 2021;Ma et al. 2017;Kapadia et al. 2020;Singh et al. 2021).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 SNR>12 as the detection criterion, where total SNR is converted from the multi-band matched filtering outputs by the analytical expression of SNR 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 multi-band localization scheme.We leave a dedicated long signal early warning pipeline and simulations with more realistic mock data to future work.

RESULTS
For each simulation, we perform multi-band matched filtering from 60 minutes before merger with low frequency cutoff at 5 Hz.Fig. 2 shows the cumulative number of events for different negative latencies.∼ 10 events can be localized within 100 deg 2 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 2 up to 6 minutes before merger.
Extreme early warnings are possible.Several events in our simulation are detected and preliminarily local- Cumulative number of detections and 90% confidence sky localization areas for the 68000 BNS simulations (roughly one month of observation).The corresponding detection efficiency is labeled in the right y axis.We choose 10 different negative latencies (from 50 minutes to post-merger) and the curves show the cumulative distribution of 90% areas of events that are detected at those times.
ized 40-50 minutes before merger and this number could be underestimated since our analysis has hard cutoffs at 5 Hz and one-hour negative latency.ET would be able to collect sensible data down to ∼ 3 Hz and trigger even earlier detections (Branchesi et al. 2023;Borhanian & Sathyaprakash 2022).However, BNS with high negative latencies are not likely to be well localized until more data comes in, bringing higher SNRs, wider frequency bands and a longer equivalent network baseline.Multi-band analysis helps update detection statistics and skymaps on-the-fly.Fig. 3 shows the evolution of skymaps and localization areas in our simulation.The example skymap 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 2 , 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 2 ∼ 20 minutes before merger in Fig. 2 are generally detected 40-50 minutes before merger.Fig. 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 multi-band localization scheme is reasonably self-consistent.The lines for 30+ minutes before merger have larger statistical fluctuations due to the insufficient number of samples.
We tested the time cost of SealGW and Bayestar calculation with the same data (full-bandwidth ET+2CE network) and skymap resulution (n side = 2048, finest 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.1s to read and preprocess the data from pipeline outputs.

CONCLUSION AND DISCUSSION
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 multi-band analysis with the fast localization algorithm For other bands, we randomly select 150 events and plot their error bar in black.The error bar is calculated from a binomial distribution, and we note that it only converges to (0%, 0) and (100%, 1) when the sample size is sufficiently large.
SealGW.We show that this is a efficient scheme for premerger localization.
Multi-band 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 2 at ∼ 10 minutes before merger.10 deg 2 can be achieved ∼ 6 minutes before merger.Since widefield optical transient facilities usually have field of view of 1-10 deg 2 (see summaries in Sachdev et al. (2020); Ronchini et al. (2022)), the precise pre-merger 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 pre-merger signals and disseminating skymaps 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 multi-band 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.
Data availability: The skymap fits files for 30, 20, 10, 6 and 1 minutes early warning and post-merger triggers are openly available in zenodo (Hu & Veitch 2023b).Further data and scripts for reproducibility are available upon reasonable request.
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.QH is supported by CSC.JV is supported by STFC grant ST/V005634/1.The matched filtering speed (detection latency) depends on the efficiency of the detection pipeline, and data conditioning usually takes ∼ 0.1s to read and pre-process the data from pipeline outputs.

Figure 1 .
Figure 1.Multi-banding scheme for this work.Left axis:An illustration of a chopped GW waveform that is alternately colored for different waveform bins and sampling frequencies.We use two-minute segments with 256 Hz sampling frequency for waveforms from 60 minutes before merger to two minutes before merger, and one-minute segments for the last two minutes with 1024 Hz and 4096 Hz sampling frequencies, respectively.Right axis: The cumulative SNR of the signal, showing the contribution of SNR from each time segment.
Figure2.Cumulative number of detections and 90% confidence sky localization areas for the 68000 BNS simulations (roughly one month of observation).The corresponding detection efficiency is labeled in the right y axis.We choose 10 different negative latencies (from 50 minutes to post-merger) and the curves show the cumulative distribution of 90% areas of events that are detected at those times.

Figure 4 .
Figure 4. P-P plot of localizations in our simulation at different bands.For 50 min (light blue, sample size = 8), 40 min (brown, sample size = 18), and 30 min (orange, sample size = 53), error bars are plotted individually with their own colors.For other bands, we randomly select 150 events and plot their error bar in black.The error bar is calculated from a binomial distribution, and we note that it only converges to (0%, 0) and (100%, 1) when the sample size is sufficiently large.

Figure 5 .
Figure 5.Time cost of running SealGW and Bayestar for ET+2CE network on a 2.44 GHz processor with different number of threads, excluding the time costs of matched filtering and data conditioning.The two algorithms are tested with the same data, and skymaps are calculated to the same level of resolution (n side = 2048, finest pixel = 0.0008 deg 2 ).The matched filtering speed (detection latency) depends on the efficiency of the detection pipeline, and data conditioning usually takes ∼ 0.1s to read and pre-process the data from pipeline outputs.
Skymap evolution.Upper panel: An example skymap for a 1.4 + 1.4M⊙ BNS at 1000 Mpc detected 30 minutes before merger with a network SNR of 12.The SNR increases to 17 at 20 minutes before merger, 31 at 10 minutes, 95 at one minute and 130 after merger.We show the 90% localization contours at different negative latencies.The injection sky location is marked with a cross.Lower panel: Evolution of 90% confidence localization areas of early warning events in our simulation.The example in upper panel is plotted in red line.