© The Authors 2023

1. Introduction

Gamma-ray bursts represent the most relativistic transient events in the sky. Discovered serendipitously (Klebesadel et al. 1973), they are now known to have an astronomical origin in the end point of massive stars, broadly classified as short duration events (SGRBs), some with extended emission (SGRBEEs), and long-duration events (LGRBs, reviewed in e.g., van Putten et al. 2014b). The end of life of relatively massive stars is believed to be the origin of neutron stars and stellar mass black holes by gravitational collapse when their nuclear burning phase ceases. Neutron stars can be found in supernova remnants harboring pulsars, as single objects such as Vela and the Crab (Hewish 1970), but also in binaries, notably the Hulse–Taylor system PSR1913+16 (Hulse & Taylor 1975) and the short-period binary PSR J0737-3039 (Burgay et al. 2003).

While pulsar emission gives unambiguous evidence of spinning neutron stars, identifying stellar mass black holes (at similar levels of rigor) appears more challenging, despite their common astronomical origin in core-collapse of massive stars. Following the seminal discovery of cosmological redshift of gamma-ray bursts in optical follow-up observations (e.g., Bloom et al. 2001) to an X-ray afterglow to GRB970228 by BeppoSAX (Costa et al. 1997), LGRBs have been identified thanks to core-collapse supernovae–energetic events of type Ib/c from relatively more massive stars (Galama et al. 1998; Hjorth et al. 2003; Stanek et al. 2003; Matheson et al. 2003; Modjaz et al. 2006; Guetta & Della Valle 2007; Kelly et al. 2008). On the other hand, GRBs originating in compact binary coalescence (Paczynski 1986) increasingly appear to be of the short-duration variety, as notably demonstrated by GW170817–GRB170817A (Abbott et al. 2017); however, this may be extended to include long-duration events, depending on black hole spin (van Putten & Levinson 2003) or progenitor binary evolution (e.g., Rueda et al. 2021).

The astronomical origin of GRBs identified with the end-point of stellar evolution leaves the central engine producing the ultra-relativistic outflows responsible for their non-thermal gamma-ray emission to be either a strongly magnetized neutron star (magnetar) or black hole (e.g., Piran 2004; Piran et al. 2019; Piro et al. 2019; Nakar 2020). In contrast to neutron stars, electromagnetic observations alone appear to be insufficient to conclusively distinguish between the two, despite a half-century of multi-wavelength observations covering GRBs by their prompt gamma-ray emission and afterglows in X-rays and optical-radio.

The recent advance of the LIGO–Virgo–KAGRA (LVK) detectors offers unprecedented and independent opportunity to advance this conundrum. Quite generally, an astrophysical black hole central engine powering an extreme transient event will be exposed to surrounding high-density matter. Gravitational radiation thereby produced is expected to be dominant over MeV-neutrinos and electromagnetic radiation when powered by it ample angular momentum energy reservoir in its angular momentum (van Putten & Levinson 2003) which, according to the Kerr metric, readily exceeds the same of a neutron star by an order of magnitude. True calorimetry hereby promises the break the degeneracy of neutron stars and black holes as the central engine, upon including total energy output in gravitational radiation (GW-calorimetry).

Observational opportunities may be found in compact binary mergers involving a neutron star (Abbott et al. 2017) and core-collapse supernovae, especially those associated with LGRBs (LSC 2018; van Putten et al. 2019b). At a distance D ≃ 40 Mpc (Coulter et al. 2017; Cantiello et al. 2018), GW170817 presents a fortuitously nearby event offering a unique prospect to witness the formation of stellar mass black holes in gravitational collapse, directly or delayed, of an intermediate hypermassive neutron star formed in the immediate aftermath of this double neutron star merger (e.g., Murguia-Berthier et al. 2021) – in addition to the new opportunities for probing physical properties of the neutron star progenitors (Abbott et al. 2018; Drago & Pagliara 2018; de Pietri et al. 2019; Bauswein 2019).

In association with GRB170817A, the precursor emission GW170817 makes it feasible to achieve a key objective of multimessenger observations (Acernese et al. 2007). An earlier demonstration of significance in timing across different observational channels is found in SN1987A, with essentially coincident arrival times of an MeV-neutrino burst detected independently by Kamiokande II and IMB, and alongside an optical light curve of SN1987A in the satellite galaxy LMC. For GW170817, independent pickup of the merger signal GW170817 by H1 and L1 is likewise coincident (within the light travel time of 10 ms between H1 and L1), followed by GRB170817A in NGC4993. At the distance of 40 Mpc, however, any further MeV-neutrino emission is a priori undetectable. Quite generally, timing coincidences carry observational significance based on a uniform prior on event timing on scales less than a Hubble time, subject only to instrumental constraints.

Gravitational wave observations of compact binary mergers is realized at observational sensitivity close to the limit defined by strain-noise amplitude of the detectors (detector limit; Abbott et al. 2021). This limit is defined by the theory of ideal matched filtering in the face of finite detector noise exploiting phase-coherence in correlation with model gravitational-wave templates over a large number of wave periods (e.g., Abbott et al. 2020). Ideally, un-modeled searches are pursued at similar sensitivities, even when phase-coherence is limited to intermediate time-scales (between the period of the wave and the total duration of the event). For transient emission featuring ascending or descending frequencies in gravitational-wave emission, broadband spectrograms provide a suitable starting point to the search for transient emission by image analysis. High sensitivity spectrograms can be realized by butterfly filtering over a dense bank of time-symmetric templates of intermediate duration τ with sensitivity exceeding that of Fourier-based spectrograms by over an order of magnitude (van Putten et al. 2014a). It and high-resolution mage analysis can be accelerated on modern graphic processor units (GPUs) with high bandwidth memory (Appendix D). This opens a window to probing the central engine of GRB170817A by GW-calorimetry at sensitivities on par with the pre-cursor GW170817.

The statistical significance of a candidate gravitational-wave event derives from the probability of false alarm (PFA), that is, Type I errors assuming the null-hypothesis of no gravitational waves passing by the LIGO detectors. Multimessenger observations involving the H1 and L1 detectors suggest exploiting both energy and event times t_k = t_k(t): a Boolean stochastic variable over {0, 1} as a function of real time t indicating an event k at some specific instant t_s (e.g., Brown et al. 2004; Donges et al. 2016) (Fig. 1). For a candidate emission feature, PFAs may derive from consistency conditions applied to a probability density (PDF) of t_s. Consistency in timing can be contextual or acontextual. Quite generally, this can be approached with and without priors from other observational channels, shown here in merged and individual H1 and L1 analyses, respectively.

Fig. 1. Schematic of event timing in the first double neutron star merger detected (top): the multimessenger event GW170817 at merger time tm with accompanying GRB170817A at tGRB = tm + tg following a gap tg ≃ 1.7 s (bottom). The central engine of GRB170817A formed in this gap by causality, whereby tm < ts < tGRB for the start-time of accompanying gravitational radiation. This gap condition poses a discrete event time problem for the probability p of false alarm (PFA) for a candidate detection to be in the 1.7 s gap between GW1709817–GRB170817A: p = tg/T marked by the (unique) global maximum of an indicator function over an observational time T given the uniform prior on astrophysical event times (Appendix A). When sufficiently large, the energy output ℰGW marked by ts can break the degeneracy between neutron stars and black holes.

Event timing of a candidate post-merger emission feature to GW170817 associated with the central engine of GRB170817A is subject to the gap condition (Fig. 1):

$$\begin{array}{c}\hfill C_1:t_sϵG,\end{array}$$(1)

given the interval G = [t_m,t_m+t_g] as an observational prior on its start-time, t_s, by causality; C₁ hereby introduces a Boolean-valued statistic, t_k(t) = 1 (t ϵ G), and 0 otherwise. Given the uniform prior on astrophysical event times, such a discrete event time carries a PFA equal to t_g/T over an observation of duration T as a conventional probability. As elucidated in Appendix A, this is analogous to the odds of winning a bet in roulette provided that the measurement of t_s is sufficiently precise.

Combining PFAs in one-class classification schemes (e.g., Simonson et al. 2017) is common in image analysis, here spectrograms in Sects. 3 and 4. Quite generally, we differentiate between continuous or discrete distributions (e.g., Block et al. 2006). The product of two p-values of continuous distributions, for instance, signal amplitude, is not a p-value, but can be merged by the celebrated Fisher’s combined probability test (Fisher 1932, 1948; Simonson et al. 2017; Lyone & Wardle 2018). For a pair of statistically independent p-values p_i (i = 1, 2), Fisher’s method defines an equivalent p-value by equivalent information content I_i = −log p_i according to ${\chi }_F^2=-2logP$, where P = p₁ × p₂ and ${\chi }_F^2$ denotes the χ² distribution with F = 4 degrees of freedom. The merged p-value derives from the cumulative density function of ${\chi }_F^2$. Fisher’s method has been widely used and indeed led to a number of potential improvements (Whitlock 2005; Heard & Rubin-Delancy 2017). When p-values p₁ and p₂ are small, p = P(1−log P) gives a merged p-value that is effectively equivalent to Fisher’s (Theiler 2004; Lyone & Wardle 2018). This is not the case in the present discrete event timing analysis of C₁, where the PFA derives as a conventional probability for Eq. (1) that can be merged as such with other PFAs (e.g., Theiler 2004).

In the next section, we zoom in on the complex merger sequence of the double neutron star merger GW170817 with the aim of identifying the central engine – a hyper-massive neutron star or rotating black hole – of the associated GRB170817A by GW-calorimetry and event timing, previously introduced in van Putten & Della Valle (2019, 2019a). These works are hereafter referred to Paper I and Paper II, respectively.

To this end, we present two statistically independent analysis of H1 and L1 data by model-agnostic spectrograms. Observational results are compared and combined based on high-resolution PDF(t_s)s realized by modern exascale heterogeneous computing (Appendix D). In particular, a joint PFA from event timing, factored over PFAs derived from PDF(t_s)s of each, provides a significantly enhanced level of confidence in the observations.

To start (Sect. 3), we revisit and extend our analysis of merged (H1,L1)-spectrograms in Paper I and II with (a) a one-parameter family of response curves calibrated by signal injection experiments and (b) PDF(t_s) generated over an extended foreground derived from a large number of small time-slides between H1 and L1, satisfying (c) a condition of extremal clustering robust against false positives from data anomalies. Statistical significance is expressed by PFA₁ from event timing by causality in the context of GW170817–GRB170817A (C₁) (Appendix A).

Next (Sect. 4), we present a statistically independent analysis of H1 and L1 individually with accompanying PDF(t_s). Statistical significance is expressed by PFA₂ of consistency in timing, t_s, by cross-correlating their respective PDF(t_s)s (C₂). In Sect. 5, we discuss our results and the observational consequences for the central engine of GRB170817A with joint PFA factored over PFA₁ and PFA₂ based on statistical independence of mean and difference in event time observations, providing a key advance over previous analysis. We summarize our conclusions (Sect. 6) with an outlook (Sect. 7) on upcoming opportunities to probe merger sequences involving neutron stars and core-collapse supernovae in the Local Universe (Sect. 6).

2. Rejuvenation in gravitational collapse

The double neutron star (DNS) merger sequence of GW170817 is complex, involving a hyper-massive neutron star (HNS) formed in the immediate aftermath (Abbott et al. 2017; Dai 2019; de Pietri et al. 2020; Murguia-Berthier et al. 2021) with the possibility of (delayed) gravitational collapse to a rotating black hole (BH; Akutsu et al. 2020):

$$\begin{array}{c}\hfill \mathrm{DNS}\to \mathrm{HNS}\to \mathrm{BH}+\mathrm{GRB}170817\mathrm{A}+\mathrm{AT}2017\mathrm{gfo}+\cdots ,\end{array}$$(2)

with the ellipses referring to additional post-merger emission in MeV-neutrinos and gravitational radiation. In a continuing gravitational collapse – direct or delayed – a black hole will form surrounded by a high-density disk or torus from debris and dynamical mass-ejecta from the merger prior (e.g., Rosswog et al. 1999; Ravi & Lasky 2014; Baiotti et al. 2008; Baiotti & Rezzolla 2017; Ciolfi 2020). In Eq. (2) the GW170817–GRB170817A association appears quite secure with independent p-values from continuous random variables: temporal (P_temp = 5 × 10⁻⁶) and spatial (P_spat = 1 × 10⁻²) consistency in the LIGO H1 and L1-detectors, on the one hand, and GRB170817A detected by Fermi Gamma-ray Burst Monitor (GBM) (Abbott et al. 2017).

In Eq. (2), the potential for breaking aforementioned degeneracy between a neutron star or black hole remnant derives, in part, from a rejuvenation of the energy reservoir in angular momentum, E_J, in the process of gravitational collapse. By conservation of specific angular momentum a/M = sin λ in prompt gravitational collapse of the HNS with spin-energy $E_J^{-}\simeq (1/2)J\mathrm{\Omega }$, $E_J^{+}=2M{c}^2{sin}^2(\lambda /4)$ of the black hole (Kerr 1963) produces:

$$\begin{array}{c}\hfill \frac{E_J^{+}}{E_J^{-}}=\frac{{(R/R_g)}^2}{10{cos}^2(\lambda /4){cos}^2(\lambda /2)}\simeq 2.5{\left(\frac{R}{18\mathrm{km}}\right)}^2.\end{array}$$(3)

Here, the right hand size refers to GW170817, using J = IΩ for the moment of inertia I ≃ (2/5)MR² of the HNS with gravitational radius R_g = GM/c², radius R and angular velocity Ω with Newton’s constant G and the velocity of light c. Further enhancement in a/M might derive from hyper-accretion (e.g., Bardeen 1970; Levinson & Globus 2013), although this does not appear to be called for in GW170817. Therefore, in the case of the gravitational collapse of a ∼3 M_⊙ hyper-massive neutron star, E_J increases by Eq. (3), setting the initial condition for potentially significant post-merger emission from a black hole-torus system.

2.1. Calorimetry in gravitational radiation

Indeed, rejuvenation Eq. (3) can safely account for Extended Emission feature in (H1,L1)-spectrogram with (Paper II):

$$\begin{array}{c}\hfill {\mathcal{E}}_{\mathrm{gw}}=(3.5\pm 1)\%M_{\odot }{c}^2.\end{array}$$(4)

Notably, (4) breaks the degeneracy between neutron stars and black holes, since ℰ_gw ≃ 6.3 × 10⁵² erg exceeds the limit of rotational energy of the hypermassive neutron star progenitor, even more so after the initial one second of neutrino cooling and internal dissipation leading to a uniform rotator – a supermassive neutron star (Beniamini & Lu 2021). Further by the observed gravitational-wave frequencies f_gw ≲ 700 Hz would imply a spin frequency well-below the Keplerian frequency limits (Haensel et al. 2009) for quadrupole emission at twice its spin-frequency, limiting its rotational energy to:

$$\begin{array}{c}\hfill E_J^{-}=\frac{{\pi }^2}{5}{f}_{\mathrm{gw}}^{2}M{R}^2\lesssim 1.6\times {10}^{52}\left(\frac{M}{2.5M_{\odot }}\right){\left(\frac{R}{18\mathrm{km}}\right)}^2\mathrm{erg}.\end{array}$$(5)

Thus, Eq. (4) exceeds Eq. (5) by a factor greater than 4, defined by maximal values of mass M and radius R. Yet, Eq. (4) is in quantitative agreement with the output of black hole spin-down against a non-axisymmetric high-density disk or torus (Paper II).

Accordingly, GW-calorimetry hints at continuing gravitational collapse to a black hole in Eq. (2). In contrast, calorimetry in the (post-merger) electromagnetic observation of GRB170817A and kilonova AT 2017gfo (Connaughton 2017; Savchenko et al. 2017; Mooley et al. 2018a,b; Ascenzi et al. 2021) shows a combined energy output limited to

$$\begin{array}{c}\hfill {\mathcal{E}}_{\mathrm{EM}}\simeq 0.5\%M_{\odot }{c}^2.\end{array}$$(6)

This output falls woefully short to break the degeneracy between a neutron star of black hole remnant in the chronicle Eq. (2), prompting the need for GW-calorimetry outlined above.

Even so, electromagnetic observations do provide us with crucial timing information. Post-merger rejuvenation in Eq. (2) appears delayed by constraints on the lifetime of the hypermassive neutron star (Pooley et al. 2018; Gill et al. 2019; Radice et al. 2018; Lucca & Sagunski 2020), inferred from a “blue” component of the associated kilonova AT 2017gfo (Smartt et al. 2017; Pian et al. 2017) (but see Ren et al. 2019; Lü et al. 2019; Piran et al. 2019 for a possibly long-lived supermassive neutron star remnant). Reviewed in Murguia-Berthier et al. (2021), this points to a lifetime of ≲1 s based on various independent studies (Granot et al. 2017; Gottlieb et al. 2018; Nakar et al. 2018; Metzger et al. 2018; Xie et al. 2018; Gill et al. 2019; Lazzati et al. 2020; Hamidani et al. 2020).

2.2. Sensitivity requirements

Performing GW calorimetry Eq. (4) to post-merger emission requires sensitivity on par with the CBC search for GW170817. Gravitational radiation from non-axisymmetric disks or tori swept up by a rotating black hole lack the phase-coherence over long duration, normally exploited in model-dependent searches applicable to compact binary coalescence (CBC).

Existing un-modeled searches by power-excess methods (Abbott et al. 2017, 2019) fall short by a threshold for detection of ℰ_th,gw ≳ 6.5 M_⊙c² (Sun & Melatos 2019). This represents a sensitivity to post-merger signals of 0.3% compared to ℰ_gW ≃ 2.5% M_⊙c² in GW170817 itself (Abbott et al. 2017), a priori in an excluded zone of the parameter space by the total progenitor mass-energy E_DNS ≲ 3 M_⊙c² of GW170817.

Here, we apply butterfly matched filtering which aims to capture signals with frequencies slowly wandering in time. Using a dense filter bank of time-symmetric templates, it realizes equal sensitivity to ascending and descending branches relevant to the present study of merger and post-merger emission to Eq. (2), squarely in the allowed zone ℰ_GW < E_DNS of any post-merger radiation.

As a pass filter for frequencies with finite time rate-of-change in frequency, constant frequency signals are relatively suppressed (van Putten 2016), while stochastic signals can be detected at sensitivity far beyond the sensitivity of conventional Fourier analysis (van Putten et al. 2014a). A precise characterization for the purpose of the present study follows and is based on detailed signal injection experiments.

3. Analysis of merged (H1,L1)-spectrograms

In this section, we consider (H1,L1)-spectrograms merged by frequency coincidences, generated by butterfly matched filtering (Papers I and II). Merged (H1,L1)-spectrograms were considered for the following steps (Appendix C): (1) the sensitivity is the same for ascending and descending chirps by time-symmetry of matched filtering templates (this step is un-modeled). (2) parameter estimation is performed by an image scan, quantifying features by a goodness of fit over a family of descending chirps using χ-image analysis (detailed in Papers I and II).

Extending previous work, we analyze merged spectrograms by application of surrogate time-slides of the following type: (S₀) foreground in conventional analysis of merged (H1,L1)-spectrograms with zero time-slide, Δt = 0; (S₁) extended foreground over small time-slides Δt less than the intermediate duration τ of our templates; (S₂) background from time-slides Δt greater than τ. We note that τ is considerably greater than the light travel time between H1 and L1.

We extract PDF(t_s) from extended foreground (S₁). This is a principle extension to previous work, restricted to foreground (S₀), that is, zero time-slide. A high-resolution PDF(t_s)s produced by heterogeneous computing (Appendix D) permit accurate parameter estimation, rigorous application of event timing (Sect. 3.2, Appendix A), robustness against false positives by extremal clustering (Sect. 3.3), and enable a direct comparison and combination of results between statistically independent analyses (Sect. 4).

3.1. Response curves to descending branches

Figure 2 shows calibrated response curves computed by model signal injections. In what follows, all injections are a model merger chirp followed by a post-merger descending chirp, where the first is fixed according to the parameters of GW170817 (Paper II, Fig. 2).

Fig. 2. Eight response curves of the two-level pipeline by χ-image analysis of merged (H1,L1)-spectrograms produced by butterfly matched filtering using a total of 56 signal injection experiments on LIGO O2 data covering GW170817 (top-left panel). These model injections are defined by Eq. (7). The dashed line is the mean of $\widehat{\chi }$-peaks averaged over (tiny) time-sides Δt within the light-travel time of 10 ms between H1 and L1 identified with extended emission to GW170817EE in a post-merger descending branch (lower panels). Precision realized is shown by parameter recovery (top-right panel). The parameter “field-of-view” covered for descending chirps covers energies 0 < ℰgw < 8% M⊙c2 (7 steps) and timescales of descent τs (eight steps, in units of seconds) with start frequency fs = 680 Hz for a model source distance D = 40 Mpc given by GW170817. These injection signals are in the shot-noise dominated regime of the LIGO detectors, where noise increases with frequency and hence sensitivity decreases with τs. Inserted is a plot of observed τs of post-merger emission to GW170817 (black square), overlaid to the calibration curve (blue curve) to infer uncertainty (along abscissa) from observed scatter (along the ordinate). Sample of a single injection of a combined merger and post-merger signal (1810–1825 s) in the snippet of H1L1-data containing GW170817 (tm = 1842.43 s) (reprinted from Paper II), shown in the lower panel. The present 56 injections Eq. (7) vary over energy ℰGW and time-scale of descent τs, keeping merger signal injection fixed at the parameters of GW170817.

Our descending chirp follows a fit to the post-merger feature attributed to radiating E_J of the compact remnant, distinct from an ascending chirp in binary coalescence, given by (Papers I and II):

$$\begin{array}{c}\hfill f(t)=(f_s-f_0)e^{-(t-t_s)/{\tau }_s}+f_0.\end{array}$$(7)

These curves introduce four parameters: a start frequency, f_s, late-time frequency, f₀, in addition to aforementioned (t_s, τ_s) used in the χ-image analysis of spectrograms. Stable estimates derive in goodness of fit over data-segments of Δt = 7 s, moving over small time-steps of 33 ms (supplementary data, Paper I).

Our response curves are generated by injection experiments cover the two-parameter space of energies and τ_s,

$$\begin{array}{c}\hfill {\mathcal{E}}_{\mathrm{GW}}=(0.5-8)\%M_{\odot }{c}^2,{\tau }_s=(1-4.5)\mathrm{s},\end{array}$$(8)

using a source distance of D = 40 Mpc of GW170817 (Cantiello et al. 2018), f_s = 680 Hz, and f₀ = 98 Hz.

We identify Eq. (7) in merged (H1,L1)-spectrograms by χ-image analysis parameterized by (t_s, τ_s, f_s, f₀) at a combined resolution of 0.5 MHz (2¹⁴ steps in (τ_s, f_s, f₀) at 30 Hz in t_s) in a high-frequency scan over 10⁹ parameter values over the 2048 s frame of H1L1-data sampled at 4096 Hz. Such is applied following a time-slides applied to H1 and L1 (prior to generating spectrograms).

Figure 2 shows results produced by 56 injections over ℰ_GW in seven steps and τ_s in eight steps. The response curves shown are averaged over time-slides covering the light-travel time of 10 ms between H1 and L1 using 41 “tiny” time-slides. (A movie of this injection analysis has been published online¹.)

In this process, an indicator function χ = χ(t_s, τ_s, f_s, f₀) represents a normalized (H1,L1)-correlation strength over strips along family of curves Eq. (7) (Paper I, supplementary data). This is reduced to $\widehat{\chi }(t_s,{\tau }_s)=\underset{(f_s,f_0)}{max}\chi (t_s,{\tau }_s,f_s,f_0)$, namely, upon projecting out (f_s, f₀) by maxima over (f_s, f₀) for each (t_s, τ_s). These projections are necessary by memory limitations, when scanning over a large number of 2 × 10⁹ parameter values per frame of 4096 s at 4096 Hz sampling rate.

Response curves vary with τ_s due to shot-noise in the H1 and L1 detectors noise above a few hundred Hz. Small τ_s corresponds to fast descent in frequency, putting injections mostly over low frequencies close to f₀, where detector noise is less and hence sensitivity is high. Conversely, large τ_s puts injections close to f_s, where shot-noise adversely affects sensitivity.

3.2. PFA from event timing based on PDF(t_s)

Significance of a candidate event – satisfying extremal clustering (Sect. 3.3, below) – is estimated by causality in formation of a central engine of GRB170917A (Papers I and II): the gap condition (1) t_m < t_s < t_m + t_g, t_g = 1.7 s, between GW170817 and GRB170817A, schematically indicated in Fig. 1.

Including causality in butterfly matched filtering with template duration τ, we consider the slightly more restrictive gap condition,

$$\begin{array}{c}\hfill {C_1}^{\prime }:\tau <t_s-t_m<t_g,\end{array}$$(9)

in the reduced gap G′ of width t_g − τ. Statistical significance is expressed by the PFA: Type I errors in $C_1^{\prime }$ according to t_s satisfying t_k(t_s) = 1.

Formally, as per Sect. 1, t_k(t):[0, T]→{0, 1} is a Boolean random variable over our snippet of H1L1-data of duration T = 2048 s. In the present case, one such candidate event (Fig. 4) marked by the (unique) global maximum at event time t = t_s in an indicator function χ(t) over [0, T], satisfying extremal clustering over extended foreground S₁. Under the null-hypothesis of noise only and no signal, t₁(t) satisfies a uniform prior on 0 ≤ t ≤ T. The Type I error, t₁(t) = 1, in finding t_s ϵ G (or G′) (satisfying C₁ (or $C_1^{\prime }$) by mere chance) is determined by G as a subset of [0, T]. Upon partitioning of [0, T] into n = 1204 intervals of size |G| = t_g, this Type I error carries a PFA given by the ordinary probability p₁ = 1/n = t_g/T. We note that this estimate is conservative compared to the same derived from G′, since |G′| = t_g − τ.

Next, we generated PDF(t_s) by application of χ-image analysis over extended foreground: the start-time t_s of a candidate emission feature according to Eq. (7) in our spectrograms. For descending chirps characteristic for spin-down of a compact object (neutron star or black hole) powering GRB170817A, Eq. (7) appears minimal in the “1+3” parameters t_s and (nuisance parameters) (τ_s, f_s, f₀). Figure 3 shows PDF(t_s) to be entirely in the t_g = 1.7 s gap between GW170817 and GRB170817A. Under the null-hypothesis – a background of uncorrelated noise in H1 and L1 – over the original snippet of T = 2048 s of H1L1-data, the gap condition C₁ is represented by a Boolean random variable t_k(t) with values that are either “true” or “false,” that is: t_k(t) = 1 if t ϵ G and t_k(t) = 0 otherwise.

Fig. 3. PDF of the start time ts (blue) of post-merger gravitational radiation to GW170817 following rejuvenation in gravitational collapse to a rotating black hole, setting the initial data to GRB170817A (top-left panel). The PDF shown (bin size 0.05 s) is superposed to estimates of the lifetime of the progenitor hypermassive neutron star formed in the immediate aftermath of the merger (reviewed in Murguia-Berthier et al. 2021). Highlighted is the 1.7 s gap (magenta) between merger time tm = 1842.43 s and the onset of GRB170817A. Post-merger emission to GW170817 appears as a descending chirp in the (H1,L1)-spectrogram merged by frequency coincidences |fH1 − fL1| < 10 Hz (lower panels). For the merged spectrogram, parameter estimation (Fig. 4) derives from time-slides less than the duration τ of the templates (11) against a background generated by larger time-slides. Post-merger gravitational radiation is identified with the formation of the central engine of GRB170817A, whose duration is constrained by the lifetime of black hole spin (bottom panel).

Given G, t_s in G according to PDF(t_s) (Appendix A) carries a PFA p₁ equal to Pr(t ϵ G) satisfying (Paper I):

$$\begin{array}{c}\hfill p_1=\frac{t_g}{T}=\frac{1.7}{2048}=8.3\times {10}^{-4}.\end{array}$$(10)

Causality over a background T hereby carries a PFA with two-sided Gaussian-equivalent significance of 3.33 σ and false alarm rate (FAR) of about 1 per month, defined by FAR = PFA/T (e.g., Williams & Huntington 2018).

3.3. Extremal clustering in extended foreground

Parameter estimation PDF(t_s) (also PDF(τ_s)) for our candidate feature is derived from merged (H1,L1)-spectrograms over an extended foreground (S₁) by application of small time-slides in the interval,

$$\begin{array}{c}\hfill I:\left|\mathrm{\Delta }t\right|<\tau ,\end{array}$$(11)

where τ ≃ 0.5 s is the intermediate duration of the templates in butterfly matched filtering. We apply this to the original snippet of T = 2048 s of H1L1-data in 99 steps symmetric about zero, including a neighborhood of zero covered by 41 tiny time-slides within 10 ms.

Following Eq. (11), the parameter estimations of candidate features give rise to “clustering” about peaks ${\widehat{\chi }}^{\ast }=\widehat{\chi }(t_s^{\ast },{\tau }_s^{\ast })$ shown in Fig. 4, where the remaining two parameters (f_s, f₀) in our χ-image analysis have been projected out by taking maxima. Such results by χ-image analysis of spectrograms for each Δt. In the extended foreground thus produced, clustering in data points – the output of χ-image analysis of merged (H1,L1)-spectrograms – is identified in the box,

$$\begin{array}{c}\hfill B:|{\tau }_s-{\tau }_s^{\ast }|<\tau ,|t_s-t_s^{\ast }|<\tau .\end{array}$$(12)

Fig. 4. Joint PDF(ts, τs) in Eq. (7) from merged (H1,L1)-spectrograms by χ-image analysis following time-slide analysis (top-left panels). Clustering is seen in the (ts, τs)-plane about a global maximum ${\widehat{\chi }}^{\ast }=8.76$ in the extended foreground over small time-slides Eq. (11). (The peak in the cluster in extended foreground (S1) is slightly larger than what is seen in conventional foreground (S0).) This cluster is extremal in the sense of Eq. (13). Tracing produces PDF’s of τs and the delay time ts − tm, where tm is the merger time of GW170817. There is notable consistency in τs with the duration $T_{90}^{8-70\mathrm{keV}}=2.9\pm 0.3$ s of GRB170817A (Pozanenko et al. 2018), expected from the timescale of black hole spin-down in interaction with surrounding high-density matter following gravitational collapse of the hyper-massive neutron star in the immediate aftermath of this merger. Clustering in τs and ts as a function of time-slide Δt in the post-merger emission to GW170817 extracted from (H1,L1)-spectrograms merged by frequency coincidences |fH1 − fL1| < 10 Hz (top-right panels). Contour plots of $\widehat{\chi }(t_s,\mathrm{\Delta }t)$ (maximal over τs, fs, f0) over all time-slides Δt ϵ [−32, 32] s (S0 − 2) in the bottom panels. The right panel zooms in on Δt ϵ [−1, 1] s, highlighting clustering within Δt ϵ [−τ, τ] set by the intermediate time scale τ ≃ 0.5 s of our template bank. Crucially, this cluster is extremal in the sense of Eq. (13) in extended foreground (S1).

For the background, the ensembles of peaks ${\widehat{\chi }}^{\ast }$ are identified over a partitioning of a data in cells of duration w. Each cell contains a total of Q = 161 data points (one for each time-slide). According to Eqs. (11) and (12), up to 99 may be clustered.

Accordingly, on the extended foreground (S₁) defined by domain (11), our two-level pipeline of (H1,L1)-spectrograms followed by χ-image analysis defines a map, whose image may be partially or wholly contained in the box (12). In case of the latter, the normalized cluster size:

$$\begin{array}{c}\hfill \eta =100\%\times \frac{\left|B\right|}{\left|I\right|},\end{array}$$(13)

defined by the ratio of counts in B and counts in I. The limiting case η = 100% will be referred to as extremal, shown in Fig. 4. (A movie of this time-slide analysis has been published online².)

Absent a signal, t_s and τ_s are uncorrelated. A uniform probability distribution of the event time t_s over the data at hand effectively prevents clustering. In response to a signal, t_s and τ_s tend to be clustered about some central value associated with a goodness of fit by Eq. (7). Such clustering is representative for correlations between H1 and L1 imparted by a gravitational-wave passing by Fig. 4. Clustering in spectrograms hereby increases with signal strength (van Putten 2017) and tends to persistent over small time-slides Eq. (11). By this property, the condition of extremal clustering provides a robust guard against false positives.

With essentially no free parameters, the extremal condition over extended foreground (S₁) provides a robust alternative to conventional thresholding in amplitude-based criteria over foreground. A common implementation of the latter is a pre-defined threshold of signal-to-noise ratio (S/N) of a continuously valued statistic.

Figure 4 with extremal clustering in (t_s, τ_s) over Eq. (11) confirms identification based on a global maximum ${\widehat{\chi }}^{\ast }=\widehat{\chi }(t_s^{\ast },{\tau }_s^{\ast })$ at Δt = 0 reported previously in Paper I, signaling Extended Emission post-merger in Fig. 3. This cluster is extremal following Eq. (13). Importantly, it is unique over the entire H1L1-data snippet of T = 2048 s following a partition over 1024 cells of 2 s cells. We hereby identify:

$$\begin{array}{c}\hfill t_s=(1843.29\pm 0.1)\mathrm{s},\end{array}$$(14)

satisfying $C_1^{\prime }$–a more precise statement of C₁. In Eq. (14) and in what follows, results are quoted by central value and STD of the PDF.

Estimating PDF(t_s) by Eq. (7) serves to produce a statistic t_s in which (τ_s, f_s, f₀) are nuisance parameters. Evaluated over aforementioned moving data-segments of ΔT obtains well-defined event timing t_s at the maximum of $\widehat{\chi }(t_s)$, similar though not identical to deriving a stable estimate of the duration of light curves of long GRBs in the face of large temporal variability by means of matched filtering (van Putten & Gupta 2009).

Given the uniform prior in timing (in the absence of astrophysical context), Eq. (9) is satisfied with associated probability of false alarm Eq. (10). Tracing the cluster of Fig. 4 generates PDF’s of t_s and τ_s, demonstrating:

$$\begin{array}{c}\hfill t_s-t_m=(0.86\pm 0.10)\mathrm{s},{\tau }_s=(2.91\pm 0.17)\mathrm{s},\end{array}$$(15)

where the first follows from Eq. (14).

Robustness of maximal $\widehat{\chi }$ conditional on extremal clustering in Fig. 4 pans out in uniqueness already in a minimally extended foreground of 7 time-slides extending over T′≃100 T in analyzing 50 frames of H1–L1 data, blindly selected to avoid biases in data-quality. This background analysis increases significance (10) in timing by causality by a factor of T′/T, that is, improving Eq. (10) by:

$$\begin{array}{c}\hfill p_1^{\prime }=\frac{t_g}{T^{\prime }}=\frac{1.7}{204800}=8.3\times {10}^{-6},\end{array}$$(16)

with a two-sided Gaussian-equivalent significance of 4.46 σ and FAR of about 1 per 782 years. This degree of robustness motivates a further analysis of timing based on individual analysis of H1 and L1 data with no regards to context. Starting point for is generating PDF’s in parameter estimation circumventing the use of time-slides, to be applied to each of the H1 and L1 data-channels.

4. Analysis of individual H1 and L1 spectrograms

Independent of context, statistical significance in timing can be evaluated according to consistency in t_s obtained from the mutually independent H1 and L1 observations:

$$\begin{array}{c}\hfill C_2:\left[t_s\right]\equiv t_{s,\mathrm{H}1}-t_{s,\mathrm{L}1}=0,\end{array}$$(17)

within the finite precision of our observation.

To evaluate C₂, we appeal to PDF₁(t_s) of H1 and PDF₂(t₂) of L1 and their cross-correlation function. The cross-correlation of the two will have a maximum, whose proximity to zero will be determined by the associated STD. As before, given a total duration of observation, T, this STD normalized to T provides us with a second PFA. This PFA from [t_s] is independent of the PFA(t_s) in the mean of H1 and L1 (Sect. 3).

Absent gravitational waves passing by, H1 and L1 output are uncorrelated. Given this null-hypothesis, mean and difference of timing (t_s,τ_s) are independent. A PDF(t_s,H1,t_s,L1,τ_s,H1,τt_s,L1) hereby factorizes as p_A(t_s,τ_s)p_B([t_s],[τ_s]). Integration over differences reduces it to p_A(t_s,τ_s); and to Eq. (10) when further integration over τ_s. On the other hand, integrating over (t_s,τ_s) reduces the same to p_B([t_s],[τ_s]) independently of the context, that is, with no regards to the multimessenger data GW170817–GRB170817.

The finite dispersion in the cross-correlation of PDF₁(t_s) and PDF₂(t_s) is expected to be quite small, albeit greater than the light travel-time between H1 and L1. In butterfly matched filtering over templates of intermediate duration τ, temporal resolution in a single cross-correlation between data and template is governed by τ. However, resolution is enhanced when averaging over a number of trials. To this end, we generated PDF’s by nm trials in a search extending over n seeds of our template banks in butterfly matched filtering, while output is branched into m spectrograms defined by a stride m. Here, template seeds refer to the master templates prior to time-slicing over τ and conversion to time-symmetric templates (van Putten et al. 2014a). As in the analysis of the previous section, these calculations are realized by modern heterogeneous computing (Appendix D).

We generate PDFs with a total of nm = 320 elements from n = 10 different template seeds and a stride m = 32. From a single template seed, the output of butterfly filtering comprises L = 𝒪(10⁸) lines of output of the form:

$$\begin{array}{c}\hfill {(i,{\rho }_1,{\rho }_2,f_{11},f_{11},f_{21},f_{22})}_k,\end{array}$$(18)

where t_i = i/4096 s is the sample time, ρ_i denotes the correlation between data and specific template for detector H1 (i = 1) and L1 (i = 2) and f_ij refers to the initial (j = 1) and final (j = 2) frequency of the template. For a stride m, this output branches into m spectrograms. Image analysis of each hereby produces a PDF of the parameters at hand, thus broadening our search by a factor nm in each of H1 and L1.

Inherent to our approach is a bank of templates that is dense rather than orthogonal, and approximate over an intermediate times scale (van Putten et al. 2014a; van Putten 2017). A sample i in Eq. (18) is hereby typically covered by a multiplicity of different templates. The average time-stepping:

$$\begin{array}{c}\hfill \mathrm{\Delta }{t}_{\mathit{ik}}=t_{i_{k+1}}-t_{i_k},\end{array}$$(19)

in t_i is correspondingly smaller than the sample time 1/4096 s of data; Δt_ik is typically zero or 1/4096 s with STD of about 40 μs. A stride m = 32 hereby covers about 7 samples. The associated STD in the mean of frequencies is similarly small, i.e., 𝒪(10⁻²) Hz. Offsets m₀ = 1, 2, …, m hereby provide exceptionally high-resolution time-stepping in our matched filtering analysis.

Figure 5 shows the results of parameter estimation in the (t_s, τ_s)-plane of H1 and L1 by χ-image analysis applied to the H1- and L1-spectrograms, generated over a large number of trials. Noticeably, our candidate signal appears stronger in H1 than in L1, indicated by $\widehat{\chi }(t_s^{\ast })$ and the peak in the H1- and L1-histograms over t_s. Furthermore, the results over T = 2048 s of H1L1-data for $\left[t_s\right]=t_{s,\mathrm{H}1}^{\ast }-t_{s,\mathrm{L}1}^{\ast }$ and $\left[{\tau }_s\right]={\tau }_{s,\mathrm{H}1}^{\ast }-{\tau }_{s,\mathrm{L}1}^{\ast }$ at local maxima of $\widehat{\chi }(t_s,{\tau }_s)$ satisfy a correlation coefficient of ρ_{[ts], [τs]} = 1.29 × 10⁻⁴, as expected from otherwise independent operation of H1 and L1. Consequently, we have the factorization p_B = p₂([t_s]) × p₃([t_s]), facilitating integrating out either [τ_s] or [t_s].

Fig. 5. Sample of ($t_s^{*},{\tau }_s^{*}$) identified in χ-image analysis of individual H1- and L1-spectrograms (top-left panels), representing two elements from a total of nm = 320 H1- and L1-spectrograms (lower panels). Peak values indicated (white rectangles) have $\widehat{\chi }=7.0$ (H1) and $\widehat{\chi }=5.3$ (L1). Over the snippet of T = 2048 s of H1L1-data, time differences $\left[t_s\right]=t_{s,\mathrm{H}1}^{\ast }-t_{s,\mathrm{L}1}^{\ast }$ and $\left[{\tau }_s\right]={\tau }_{s,\mathrm{H}1}^{\ast }-{\tau }_{s,\mathrm{L}1}^{\ast }$ identified by local maxima of $\widehat{\chi }(t_s,{\tau }_s)$ over the 64 segments of 32 s and 320 trials of H1- and L1-spectrograms are, as expected, uncorrelated, evidenced by a correlation coefficient ρ[ts], [τs] = 1.29 × 10−4. (Right-top-panels) PDF1(ts) and PDF2(ts) of H1 and, respectively, L1 are centered around ts − tm = (0.92±0.09) s and τs = (3.00±0.09) s. Binning over δtres = 0.025 s, cross-correlation of the two PDFs (over each segment covering H1L1-data over T = 2048 s) shows a global maximum in segment 58 containing GW170817 consistent with zero time difference within δtres = 0.025 s.

Figure 5 shows PDF₁(t_s) and PDF₂(t_s) thus extracted from H1 and L1 with their cross-correlation obtained by time-slide analysis applied to this pair. A key finding is that these PDF’s are maximally correlated at a time-slide consistent with zero within the time-slide resolution δt_res = 0.025 s. This maximum is global (maximal over all such maxima identified in the 64 segments of 32 s covering our analysis of the H1L1 data of duration T = 2048 s). It thereby carries a PFA:

$$\begin{array}{c}\hfill p_2=\frac{4\delta t_{\mathrm{res}}}{T}=4.88\times {10}^{-5},\end{array}$$(20)

with equivalent two-sided Gaussian-equivalent significance of 4.06 σ and FAR of ∼1 per 1.4 year.

In deriving Eq. (20) from |[t_s]| < δt_res at the peak of cross-correlating PDF(t_s,H1) and PDF(t_s,L1), a factor of 2 derives from absolute value of time-differences [t_s] = |t_H1 − t_L1| in H1 and L1; another factor of 2 derives from the top-hat distribution in [t_s] based on uniform priors in t_H1 and t_L1. Similar to p₁, consistency in the PDF’s of t_s is hereby formally cast under the null-hypothesis by the Boolean random variable Y ϵ{0, 1} with Y = 1 if |[t_s]| < δt_res, Y = 0 otherwise.

The mean, now from the pair of PDFs of t_s from the independent H1- and L1-analyses, shows a delay time in gravitational collapse

$$\begin{array}{c}\hfill t_s-t_m=(0.92\pm 0.09)\mathrm{s},\end{array}$$(21)

consistent with Eq. (15) based on the mean t_s in merged (H1,L1)-spectrograms.

In astrophysical context, Eq. (21) is consistent with the lifetime $0.{98}_{-0.26}^{+0.31}$ s of the progenitor hypermassive neutron star derived from jet propagation times and the mass of blue ejecta (Gill et al. 2019) and its survival time, based on neutrino cooling and internal dissipation of rotational energy (e.g., Beniamini & Lu 2021). For τ_s, the PDFs associated with the cluster (21) show τ_s = (3.1±0.1) s (H1) and τ_s = (2.9±0.2) s (L1) with mean:

$$\begin{array}{c}\hfill {\tau }_s=(3.00\pm 0.09)\mathrm{s},\end{array}$$(22)

consistent with Eq. (15) in the previous analysis of merged (H1,L1)-spectrograms. Identified with the lifetime of black hole spin, it appears to constrain the duration $T_{90}^{8-70\mathrm{keV}}$ of GRB170817A (Pozanenko et al. 2018). Thus, Eqs. (21) and (22) both appear to possess natural astrophysical context in the electromagnetic spectrum.

5. Discussion

The complex merger sequence Eq. (2) is probed by multimessenger calorimetry and event timing (Fig. 1). Following calibration by signal injection experiments (Fig. 2), a post-merger descending chirp characteristic for spin-down of a compact remnant is analyzed for output ℰ_GW and statistical significance by event timing subject to two consistency conditions, C_1, 2, on time-of-onset t_s seen in merged and individual H1- and L1-spectrograms (Sects. 3 and 4). Figure 6 shows a detailed summary.

Fig. 6. Chronicle (2) showing detailed event timing in Fig. 1. Highlighted is ℰGW by rejuvenation in gravitational collapse of the hyper-massive neutron star formed in the immediate aftermath of the merger at ≲ts in the 1.7 s gap G between GW170817–GRB170817A, in a Kerr black hole with spin-energy $E_J^{+}$ exceeding $E_J^{-}$ in its progenitor according to (3). $E_J^{+}$ amply accounts for ℰGW = (3.5 ± 1)% M⊙c2, breaking the degeneracy between neutron stars and black holes.

ℰ_GW in Eq. (4) suffices to break the degeneracy between a hyper-massive neutron star and a black hole by exceeding limits in E_J of the first. Yet, it is amply accounted for by E_J following (delayed) gravitational collapse to a rotating black hole in the aftermath of GW170817 (Sect. 2). Statistical significance is expressed by a PFA factorized over two statistically independent PFAs Eq. (10) in Sect. 3 and Eq. (20) in Sect. 4, derived from consistency conditions on effectively the mean (C₁) and, respectively, a difference (C₂) in start-time, t_s.

5.1. Insensitivity to numerical parameters

Event timing analysis is rather insensitive to numerical parameters given that certain discrete conditions are met in PDF(t_s) when evaluating C₁ (Appendix A) and in considering the difference [t_s] in C₂ as follows. PFA p₁ from C₁ (more conservative than PFA from $C_1^{\prime }$) is determined by G relative to the observational window [0,T], given that PDF(t_s) is well within G (Fig. 3, Appendix A). p₁ is hereby insensitive to trial factors that might otherwise be encountered when inferred from a statistic in continuous stochastic variables such as the S/N. It likewise does not involve Occam factors, common to likelihood analysis of the same in Bayesian approaches. The PFA p₂ from C₂ is derived from cross-correlating PDF(t_s) based on two individual H1 and L1 analyses. Sampled by nm = 320 trials, they appear to be close to the noise limit of the data. Moreover, any systematic uncertainty in PDF(t_s) cancels out in the difference [t_s] (Sect. 4).

An additional safeguard derives from extremal clustering (100% in Eq. 13). It implies the width of PDF(t_s) to be much smaller than the gap size, t_g, of GW170817–GRB170817A (Figs. 3 and Eq. (A.2)). Over extended foreground, it derives as a goodness of fit over (t_s, τ_s, f_s, f₀), representing an essentially minimal number of “1+3” parameters for descending chirps. (By count this is the same in a canonical CBC search, that is, taken over two masses and orbital ellipticity for merger times 0 ≤ t_m ≤ T.) There is no fine-tuning in PDF(t_s) derived from Eq. (7) with maxima $\widehat{\chi }=\widehat{\chi }(t_s)$ upon projecting out the nuisance parameters (τ_s, f_s, f₀).

Our focus on Eq. (2) gives a concrete example of multimessenger calorimetry and event timing analysis, summarized in the chronicle Fig. 6. The joint PFA p₁p₂ = 4.1 × 10⁻⁸ from (discrete) event timing is noticeably on par with the PFA 8.91 × 10⁻⁷ of GW170817–GRB170817A, merging the two p-values of (continuous stochastic variables) in temporal and spatial agreement between H1, L1, and the Fermi GBM reported by LIGO–Virgo (Abbott et al. 2017). A further third PFA derives from consistency in τ_s from the individual H1 and L1 observations (Table 1). In this regard, our total PFA by t_s is conservative.

5.2. Advances over previous results

We present two independent analysis of merged and individual H1 and L1 spectrograms, results of which may be compared for consistency and combined to enhance confidence levels on the basis of high-resolution PDFs (Sects. 3 and 4). To be precise, in the absence of a gravitational-wave signal, data from H1 and L1 are statistically independent (the working hypothesis in multi-detector GW analysis). Under the null-hypothesis, merged and individually, H1 and L1 satisfy uniform priors on event timing. We consider the mean (Sect. 3) and difference (Sect. 4) in event timing expressed in a statistic t_s, representing two linearly independent combinations of the H1 and L1 time-series equivalent to a rotation over π/4 in the plane. By unitary of rotations, statistical independence of the H1 and L1 data is preserved in the two statistics t_s of Sects. 3 and 4. Table 1 lists the measurement results.

A comparison with previous results derived from foreground (S₀, time-slide zero) is also opportune. Our PDF(t_s) produces t_s ≃ 0.92 s (Table 1), revising t_s ≃ 0.67 s (Paper II).

This revised estimate is illustrative for a limitation of foreground in representing a single sample in PDF(t_s) derived from extended. In the face of finite scatter seen in the plot of t_s versus time-slide Δt (Fig. 4, top right panel), foreground estimates incur a statistical uncertainty on the order of the width of PDF(t_s). Additionally, foreground incurs a systematic error given by the finite difference in signal arrival time between the two detectors.

The PDF(t_s)s from merged and individual H1 an L1 analysis permit a direct comparison between the two for their implied PFAs p₁ and p₂ (Table 1). The inferred PFAs are qualitatively distinct, however, being contextual (by the gap t_g = 1.7 s between GW170817–GRB170817) and, respectively, acontextual. Strictly speaking, a direct numerical comparison of the two is not opportune.

Nevertheless, our p₁ > p₂ may appear paradoxical in light of a visibly higher contrast of the descending branch in merged rather than individual spectrograms. Indeed, the width in PDF(t_s) is relatively smaller, approximately by a factor of $\sqrt{2}$ (Table 1), as expected. Increasing the time of observation by ×100, $p_1^{\prime }<p_2$ (Eq. (16)) with confidence level 4.46σ. In fact, including amplitude information, p₁ decreases below p₂ to 2.5 × 10⁻⁵ (4.2σ, supplementary data to Paper I). The same applied to $p_1^{\prime }$ provides a combined confidence level of 5.16σ, but this avenue is not pursued further here.

Discrete event timing producing PFAs differently from amplitude-based PFA (Appendix A). Figure 1 shows the candidate signal to be sufficiently strong for PDF(t_s) to be well within the gap of GW170817–GRB170817A, satisfying extremal clustering (Fig. 4). It hereby satisfies causality–a Boolean valued conclusion. True, in the case at hand, fixes PFA₁ at p₁ = t_g/T, where T is the time of observation. A further increase in signal strength, narrowing PDF(t_s), will not change this value.

For uniformity of presentation, PFAs are here given by event timing only, sufficient to derive a significantly improved joint PFA, factored over the independent PFA₁ and PFA₂, to validate the descending branch.

5.3. Black hole spin-down

We identified our estimate of τ_s ≃ 3 s in Eqs. (15) and (22) with the Kelvin-Helmholtz time-scale of black hole spin-down against high-density matter through a torus magnetosphere (van Putten & Levinson 2003; van Putten 2015a):

$$\begin{array}{c}\hfill {\tau }_H=\frac{E_J}{\left|{\dot{E}}_J\right|}\simeq 3\mathrm{s}{\left(\frac{\sigma }{0.05}\right)}^{-1}{\left(\frac{z}{6}\right)}^4\left(\frac{M}{3M_{\odot }}\right),\end{array}$$(23)

where σ = M_T/M in the right hand-side refers to the mass-ratio of torus-to-black hole with radius zR_g, R_g = GM/c², specialized to the present GRB170817A. This timescale Eq. (23) is derived from the first law of thermodynamics.

For a black hole with mass, M, angular momentum, J_H, and angular velocity, Ω_H, interacting with a torus magnetosphere rotating at the angular velocity Ω_T of the inner face of a surrounding torus, τ_H in (23) largely represents the dissipation, D_H of E_J, on the event horizon, leaving a net luminosity L_H = −Ṁ ≪ D_H (van Putten & Levinson 2003; van Putten 2015a):

$$\begin{array}{c}\hfill L_H=-{\mathrm{\Omega }}_H{\dot{J}}_H-D_H=L_j+L_T,\end{array}$$(24)

in the limit of small Reynolds stresses in suspended accretion, where L_j ≪ L_T refers to the luminosity L_j in a baryon-poor jet (BPJ) along an open magnetic flux tube subtended over a finite polar angle θ_H ≪ π/2 on the event horizon of the black hole and a luminosity $L_T\simeq -{\mathrm{\Omega }}_T\dot{J}$ in multimessenger radiation from the surrounding inner disk or torus. Here, the inequality L_j ≪ L_T ≃ L_GW is seen in discrepant energies in electromagnetic Eq. (6) and, respectively, gravitational radiation Eq. (4).

GRB170817A derived from a BPJ is included in Eq. (6). While GRB170817A is negligible in the total energy budget $E_J^{+}$ (Sect. 2), it provides crucial timing information t_GRB and $T_{90}^{(8-70)\mathrm{keV}}$. It can be attributed to L_j derived from a Faraday-induced potential (van Putten 2000; van Putten & Levinson 2003),

$$\begin{array}{c}\hfill U=\omega J_p,\end{array}$$(25)

of charged particles with angular momentum J_p = eA_φ along surfaces of constant flux A_φ = Φ/2π subject to frame-dragging by a black hole in its lowest energy state. While L_j for θ_H ≪ π/2 is relatively small, it may have observational consequences for GRBs and UHECRs alike and especially so when intermittent (van Putten & Gupta 2009; Shahmoradi & Nemiroff 2015; van Putten 2015b; Gottlieb et al. 2021). For a discussion on high-energy emission from black holes that are not in the lowest energy state, we refer to Rueda et al. (2022).

6. Conclusions

From analysis of H1 and L1 observations in merged and, new in this work, individual detector spectrograms producing statistically independent PFAs p₁ and, respectively, p₂, from their respective high-resolution PDF(t_s)s of an extended emission feature, a number of findings emerge:

1. A black hole central engine to GRB17017A based on ℰ_GW ≃ 3.5% M_⊙c² in a descending chirp exceeding the limits of the HNS in the immediate aftermath of GW170817 (Eqs. (4), (7), (5)), attributed to rejuvenation of E_J in gravitational collapse (Sect. 2);

2. Statistical significance with combined PFA of 4 × 10⁻⁸ (equivalent Z-score 5.48) factored over two independent PFAs in Eqs. (10) and (20) derived from event timing subject to C_1, 2 (Table 1);

3. A delay time t_s − t_m ≲ 1 s in post-merger gravitational collapse in agreement with independent estimates of the lifetime of the hypermassive neutron star (Fig. 3, Table 2), ruling out a BH-NS merger progenitor (Coughlin & Dietrich 2019);

4. A timescale of descent τ_s ≃ 3 s identified with the lifetime of black hole spin in agreement with the duration $T_{90}^{8-70\mathrm{keV}}=2.9\pm 0.3$ s of GRB170817A (Pozanenko et al. 2018). We refer in particular to Figs. 4, 6, and Table 2;

5. The gravitational-wave energy emitted in the combined ascending-descending chirp of GW170817 is about 2% of the total mass-energy of the system.

7. Outlook

Planned LVK observational runs O4-5 offer potentially important new observational opportunities to probe merger sequences involving neutron stars and central engines of energetic core-collapse supernovae in the Local Universe. EM-GW observations involving GW-calorimery and event timing offer some new tools to identify their nature (Cutler & Thorne 2002), particularly when breaking the degeneracy between (super- or hyper-) massive neutron stars and black holes.

Since black holes have no memory of their progenitor except for total mass and angular momentum, such appears notably opportune for type Ib/c supernovae as the parent population of normal long GRBs (van Putten et al. 2019b) and superluminous supernovae (Dong et al. 2015). While rare, failed GRB-supernovae nevertheless may be luminous in gravitational radiation and more frequent than double neutron star mergers.

For the planned LVK runs O4-5, core-collapse supernovae may be probed in blind or optically triggered searches in the Local Universe over distances comparable to GW170817. If detected, gravitational-wave emission is likely to identify their enigmatic central engine, significantly complementing our understanding of core-collapse events currently limited to SN1987A. In the more distant future, searches for broadband gravitational-wave radiation may put rigorous and model-independent bounds on (non-axisymmetric) mass-motion around supermassive black holes such as SgrA* by the planned Laser Interferometric Space Antenna (LISA; van Putten et al. 2019b).

Acknowledgments

The authors gratefully acknowledge a detailed reading and constructive comments from the anonymous reviewer and M.A. Abchouyeh, which greatly contributed to clarity of presentation. The first author gratefully acknowledges stimulating discussions with Gerard’t Hooft over a Nico van Kampen Colloquium on GW170817 at ITP, University of Utrecht (2019). LIGO O2 data are from the LIGO Open Science Center of the LIGO Laboratory and LIGO Scientific Collaboration (LSC), M. Vallisneri et al., 2014, Proc. 10th LISA Symp., University of Florida, Gainesville (May 18–23), [arXiv:https://arxiv.org/abs/1410.4839], funded by the U.S. National Science Foundation. The original GW170817 2048 s data set is 10.7935/K5B8566F of the LIGO Laboratory and the LSC. Virgo is funded by the “French Centre National de Recherche Scientifique (CNRS)”, the “Italian Instituto Nazionale della Fisica Nucleare (INFN)” and the “Dutch Nikhef”, with contributions by Polish and Hungarian institutes. Additional support is acknowledged from MEXT, JSPS Leading-edge Research Infrastructure Program, JSPS Grant-in-Aid for Specially Promoted Research 26000005, MEXT Grant-in-Aid for Scientific Research on Innovative Areas 24103005, JSPS Core-to-Core Program, Advanced Research Networks, and the joint research program of the Institute for Cosmic Ray Research. Computations have been performed on a dedicated platform by synaptic parallel computing for dynamical load balancing. This research is supported, in part, by NRF of Korea Nos. 2018044640 and 2021K1A3A1A16096820. M.D.V. acknowledges support from PRIN-MIUR 2017 No. 20179ZF5KS. The data underlying this article were accessed from the LIGO Open Science Center, https://www.gw-openscience.org/about/, specifically the 2048 s data set 10.7935/K5B8566F containing the merger GW170817. Derived data generated in this research will be shared on reasonable request to the corresponding author.

References

Appendix A: Discrete event timing versus S/N

GW170817-GRB170817A provides a multimessenger context by the gap G of t_g = 1.7 s in between these two events with merger time t_m and, respectively, the time-of-onset t_GRB (Fig. 1).

For a candidate post-merger emission feature associated with the central engine of GRB170817A, a PFA derives from event timing as an alternative to conventional signal-to-noise (SNR) analysis under the null-hypothesis of stationary detector noise and the uniform prior of astrophysical event timing (H0).

To illustrate this, consider an indicator function χ(t) marking the event time t_e of a candidate feature by its global maximum over a finite duration of observation T. Since χ(t_e) is continuous stochastic observable, t_e is well-defined and unique. Thus, t_e tends to be uniformly distributed over [0, T], namely:

$$\begin{array}{c}\hfill {\mathrm{Pr}}_0=p(t_e|H0)=\frac{1}{n},\end{array}$$(A.1)

with event time now discrete over n = T/t_g bins. The gap condition (1) hereby carries a PFA equal p₁ = t_g/T.

In contrast, gravitational-wave emission from the putative central engine of GRB170817A (H1) changes our expectation to Pr = p(t_e|H1) with a bias to satisfying Eq. 1, provided χ(t) is suitably devised to measure the start-time t_e = t_s of this emission. As may be expected, p(t_e|H1) > p(t_e|H0) for t_e ϵ G depends on the S/N.

Figure A.1 shows an elementary Monte Carlo simulation of the probability Pr(t_e ϵ G) of satisfying the gap condition Eq. 1 as a function of S/N, based on event times, t_e, defined by maxima of χ(t) given by the sum of normally distributed noise over n = 37 bins plus a signal at bin 30.³ For n = 37, this PFA equals the chance of winning a single bet in European roulette conform the probability theory of Blaise Pascal.

Fig. A.1. MC simulation on discrete event timing over n = 37 cells associated with a gap condition at cell 30, shown in magenta (top and middle panels). The event time, te, of a single run is the location of the maximum of an indicator function, χ(t), over all cells. A signal from a source in the gap (H1) increases the probability, p(te|H1), of detection at this location, above the PFA p(te|H0) = 1/n of stationary background noise (H0). The probability p(te|H1) increases monotonically with S/N, here derived from 1024 runs (lower panel).

A PFA from Eq. 1 rather than based on the S/N has the advantage of being independent of trial factors, though the result is a priori limited by G (relative to T) and the implementation requires PDF(t_e) to be sufficiently narrow relative to G, that is:

$$\begin{array}{c}\hfill {\sigma }_{t_e}\ll t_g\ll T.\end{array}$$(A.2)

Once (A.2) is satisfied, p₁ = Pr₀ is independent of the width of PDF(t_s), which might otherwise depend on trial factors. In this sense, (A.2) represents a discrete timing condition. As an ordinary probability, furthermore, p₁ readily combines with other PFAs.

Appendix B: Whitening H1L1-data

As a pre-processing step applied to LIGO strain data, we applied a band pass filter of 10-1700 Hz followed by whitening to suppress numerous lines, mostly violin modes which appear prominently in the spectrum of the detectors (Fig. B.1) (Abbott et al. 2020).

Fig. B.1. (H1L1-detector spectral energy density Sn(f) for the T = 2048 s data snippet containing GW170817 following a band pass filter 10-1700 Hz (left panel). Violin modes appear prominently associated with the suspension of optics. Frequencies up to about 1700 Hz can be used in injection experiments. H1 and L1 detector noise is very similar though not identical during the GW170817 event. Whitening of H1 and L1 strain data by normalisation (B.1) produces flat spectra with lines suppressed (right panels).

Figure B.1 shows the computation by the standard Welch method. Here, the frequency resolution is slightly less than in Fig. 1 of Paper II by a different partitioning in the time domain. By the time-frequency uncertainty relation, lines hereby vary in width yet with the same total energy by Parseval’s Theorem. We note, however, that constant frequencies are suppressed in our butterfly matched filtering (Appendix C).

To this end, strain data are normalized in Fourier domain by amplitude following a partitioning of the spectrum over intervals of B [Hz] (supplementary data, Paper I). Specific to the T = 2048 s snippet of H1L1-data covering GW170817 sampled at a n_s = 4096 Hz, we have N = 2¹²[T/s] = 2²³ Fourier coefficients c_k (k = 1,2,⋯N). Its Fourier spectrum is partitioned up to the Nyquist frequency $\frac{1}{2}{n}_s$ by intervals B_s: (s − 1)M < k ≤ s × M of size B[Hz], each comprising M = (B/n_s)N = BT coefficients; to exemplify this, M = 4096 for B = 2 Hz. Normalization to an essentially flat spectrum is obtained by:

$$\begin{array}{c}\hfill C_k=\frac{c_k}{A_s}\left(kϵB_s\right),\end{array}$$(B.1)

over all intervals $s=1,2,\cdots ,\frac{1}{2}{n}_s/B$, where A_s = M⁻¹∑_kϵBs|c_k| is the mean of the absolute values of c_k in B_s with |B_s| = B. Crucially, (B.1) preserves phase in each of the Fourier coefficients. Effectively the same whitening obtains upon setting A_s equal to the standard deviation of the c_k in B_s.

This procedure suppresses violin modes in the LIGO detectors, provided B is about 1-10 Hz - B must exceed the maximal width of the violin modes yet be relatively modest to effectively suppress the same. The result is essentially flat spectra (Fig. B.1), while preserving signals of interest such as the relatively long-duration merger signal GW170817, evidenced in Fig. 3. Whitening (B.1) is effective also in rendering GW170817 to be directly audible, for instance, using the MatLab function sound.

In the present application to the H1L1-data covering GW170817, we empirically verified in previous work (Paper I) that all results remain essentially unchanged for B = 2 − 16 Hz.

Appendix C: Butterfly matched filtering

Butterfly filtering is an essentially linear filter for signals with frequencies slowly wandering in time parameterized by δ (Fig. C.1). This is devised by matched filtering over an effectively dense bank of time-symmetric chirp-like templates of duration τ, intermediate to the period and the total duration of a candidate signal (van Putten 2017).

Fig. C.1. Time-frequency domain illustrating filtering signals with frequencies slowly wandering in time, that appear as trajectories with finite time-derivative (top panel). Schematically, these trajectories passing through “butterflies" (green) defined by a finite slew rate (time rate-of-change) in frequency |df(t)/dt| ≥ δ > 0 (insert). It is realized by matched filtering over a dense bank of chirp-like time-symmetric templates of intermediate duration. (Reprinted from van Putten (2016).) Identification of broadband Kolmogorov spectrum in an ensemble averaged spectrum to the Nyquest frequency of 1 kHz of 42 bright long GRBs in the BeppoSAX catalog by butterfly filtering over a dense bank of 8.64 million templates, shown as a purple line (bottom panel). This result demonstrates sensitivity to turbulence better than an order of magnitude compared to conventional Fourier analysis (blue). Note: reprinted from van Putten et al. (2014a).

In the time-frequency diagram, the passing of such signals can be schematically indicated by "butterflies" (Fig. C.1), indicated by a finite slew rate (time rate-of-change) in frequency:

$$\begin{array}{c}\hfill \left|\frac{df(t)}{\mathit{dt}}\right|\ge \delta ,\end{array}$$(C.1)

for some choice of δ > 0. Constant frequency signals fail to pass through the same, indicating relative suppression. Butterfly matched filtering is hereby distinct from Fourier analysis, that favors passing signals with relatively constant frequency.

Butterfly filtering was originally developed to extract broadband spectra of light curves of long GRBs of the BeppoSAX catalog (van Putten et al. 2014a). It identifies Kolmogorov spectrum which extends up to the Nyquist frequency of 1024 Hz of the BeppoSAX sampling rate of 2048 Hz (during the first 8 seconds of long bursts), demonstrating a sensitivity one order of magnitude beyond what is attained by conventional Fourier analysis (Fig. C.1).

The template bank of butterfly matched filtering densely covers a two-parameter range in frequency and time rate-of-change of frequency. These templates are produced by time-slicing of a seed template (Fig. C.2) over aforementioned duration τ. Time-symmetric templates obtain by linear combination with their time-reverse, realizing equal sensitivity to signals whose frequency increases or decreases in time.

Fig. C.2. Master template serving as a seed to the construction of a template bank, given by an exponential decay in frequency f0 ≥ f(t) > f1 over 0 ≤ t ≤ T (top panel). Time-slicing produces chirps of intermediate duration P ≪ τ ≪ T, where P denotes the characteristic period of the (candidate) signal of interest. Time-symmetric templates obtain by linear superposition with their time-reverse (right), realizing equal sensitivity for ascending and descending chirp-like signals. (Reprinted from van Putten et al. (2014a).) Illustration of the suppression of constant frequency input (lower panel) in the form of the application of butterfly filtering to unwhitened LIGO data (left). Lines such as shown in Fig. B.1 fail to correlate with butterfly templates, leaving a smooth and approximately white background noise (right). Note: reprinted from van Putten (2016).

Butterfly filtering is ported to LIGO data analysis with no principle change in algorithm. An early demonstration shows the suppression of lines in unwhitened LIGO data, as constant frequency signals fail to correlate with chirp-like templates (Fig. C.2).

Analysis of large data-sets using a dense bank typically containing 2¹⁸ templates requires heterogeneous computing, off-loading matched filtering to graphics processor units (GPUs) with additional procedures to manage limitations of bandwidth between GPU and host CPU given exascale computations involved (van Putten 2017, Appendix C).

Appendix D: Heterogeneous Parallel Computing

The following gives a brief description of the pipeline to probe the merger sequence Eq. 2 optimized for high throughput computation (HTC) to facilitate the analysis in §3-4 by modern heterogeneous computing. Table D.1 lists some of the relevant symbols.

D.1. Exascale compute requirements

To illustrate computational requirements, consider the snippet of 2048 s of H1L1-data partitioned over 64 segments of 32 s duration. A single spectrogram is produced by 64 × 2¹⁹ ≃ 3.4 × 10⁷ (H1,L1)-correlations over a bank with 2¹⁹ = 524288 templates. Evaluated by FFT with n = 2¹⁷ samples per segment at 4096 Hz, this comprises about 373 teraFLOPs (floating point operations) in single precision (SP, f32). Extended over 161 time-slides Δt in S_0 − 2 (§3) is aggregated to 60 petaFLOPs.

The response curves in Fig. 2 are extracted from a data-base of 56 injection experiments produced in about 56 × 60 ≃ 3.3 exaFLOPs. A similar computational effort holds for χ-image analysis applied to relatively dense scatter plots of spectrograms.

An analysis was carried out on a heterogeneous platform of multi-GPU nodes with synaptic parallel processing (below) over multi-LANs for dynamical load balancing, realizing the required high throughput computations (HTC).

D.2. Computing steps

The two-level pipeline is implemented by heterogeneous computing:

1. Pre-processing: H1 and L1 data-frames from the LIGO Open Science Center (LOSC) are paired by GPS start-time, whitened and stored as complex numbers in (H1,L1)-data in single precision in Fortran binary format (little-endian).

2. Whitening: applied to remove numerous violin modes in the LIGO data. This can be effectively carried out in the Fourier domain by normalizing the spectrum of strain data over intervals of relatively modest bandwidth B (Appendix B, Paper I);

3. Spectrograms: generated on the basis of (H1,L1) data via butterfly matched filtering is implemented in OpenCL (clButterfly) on multi-GPU nodes with high bandwidth memory (HBM) for fast evaluation in the Fourier domain, using of pre- and post-callback functions to circumvent PCIe bandwidth (van Putten 2017);

4. Rendering: applied to plot spectrograms in MatLab;

5. Parameter estimation: χ-image analysis⁴ of both merged and individal H1 and L1-spectrograms is implemented in OpenCL (clChi).

6. Post-processing: gathering results into database for statistical analysis is carried out on CPUs in MatLab.

In §3, the above is iterated over time-slides (S_0 − 2). In §4, iteration is over n template seeds and m offsets by embarrassingly parallel computing.

D.3. Limits on performance

Speed-up by off-loading tasks to GPUs depends crucially on arithmetic intensity. Matched filtering by correlations carried out in the Fourier domain tend to be low arithmetic intensity for arrays of relatively large size, when Global Memory calls are required for matrix transpose. Performance is then well below theoretical peak performance effectively limited by memory bandwidth (Fig. D.1), pointing to the need for high bandwidth memory (HBM) and a platform comprising multiple multi-GPU nodes. Further optimization obtains by circumventing bandwidth limitations of the peripheral computer interface express (PCIe) with pre- and post-callback functions in the OpenCL implementation clFFT (van Putten 2017).

Fig. D.1. Performance overview of GPU-accelerated butterfly matched filtering on a quad-GPU node. Performance of clFFT in CSP/Out-of-Place as a function of array size N for two batch sizes measured by allocation M in global memory (left panel). Array sizes exceeding 212 cause a drop in performance as matrix transpose exceeds 32 kB of local memory, limited now by 1TB/s bandwidth to Global Memory in HBM2. We note that clButterfly uses a default size N = 217 (32 s of data at 4096Hz sampling rate) with batch size over n = 128 per frame of 4096 s. A quad-GPU node hereby features a maximum of about 5 teraFLOPs, circumventing limitations of PCIe bandwidth optimized by pre- and post-callback functions (van Putten 2017). The throughput of this quad-GPU node is shown up to six times real-time over a template bank normalized to size 220 (right panel). GPU-CPU communication passes through a post-callback function, limiting data-transfer to tails exceeding a threshold κ = 2 in normalized butterfly matched filtering output, effectively 0.1% of all matched filtering results. Throughput peaks in multi-frame analysis at a cross-correlation rate of about 200kHz.

In contrast, clChi can be implemented in Local Memory, effectively realizing the computing limit performance on the same GPUs. We note, however, that scanning over (τ_s, f_s, f₀) for each step in t_s (about 2 × 10⁹ parameter values per data-frame of 4096 s) represents multiple instruction single data (MISD) processing, rather than SIMD in conventional graphics processing. Such requires care to avoid overheating the memory bus to Global Memory, by (pseudo-)randomization over t_s in distributing work Items over the large number of stream processors in the GPU. Programs clButterfly and clChi are written in C++/F90 with kernels in C99.

D.4. Multi-GPU node benchmarks

Benchmarks of butterfly matched filtering can be normalized to clFFT performance⁵ as follows. With correlations between data and templates evaluated in the Fourier domain, throughput critically depends on the efficiency,

$$\begin{array}{c}\hfill {\eta }_0=\frac{\mathrm{cross}-\mathrm{correlation}\mathrm{rate}}{\mathrm{FFT}\mathrm{transform}\mathrm{rate}}.\end{array}$$(D.1)

Due to limited bandwidth of the Peripheral Computer Interface express (PCIe) between GPU and CPU, there are variations with regard to parallel processing over multi-frame jobs (van Putten 2017). It does, however, tend to improve with larger job sizes. Net throughput can be expressed as

$$\begin{array}{c}\hfill {\mathrm{throughput}}^{-1}=\frac{\mathrm{wall}\mathrm{clock}\mathrm{time}[\mathrm{s}]}{\mathrm{data}\mathrm{length}[\mathrm{s}]\times \mathrm{template}\mathrm{bank}\mathrm{size}[\mathrm{M}]}.\end{array}$$(D.2)

Here, the template bank size is expressed in units of 1M templates. A throughput in excess of unity defines faster than real-time analysis.

Fig. D.1 illustrates D.2-D.1 for a quad-GPU configuration of an air-cooled node with AMD Radeon VII GPUs featuring η ≃ 47% and a cross-correlation rate up to 200 kHz between data and templates over segments of 32 s.

D.5. Load balancing by synaptic processing

Embarrassingly parallel computing for tasks comprising clButterfly (step 3) or clChi (step 5) in §3-4 are specified as a line-items for processing on a multi-LAN heterogeneous computing platform.

We used synaptic parallel processing in which line-items are requested by the GPU-nodes, rather than issued by a server. Before execution, these line-items are updated with the URLs for input and output data on the LAN hosting a node. This process ensures maximal throughput on heterogeneous platforms of nodes with different performance characteristics. Synaptic parallel processing is written in bash for Linux and OSX.

All Tables

All Figures

Fig. 1. Schematic of event timing in the first double neutron star merger detected (top): the multimessenger event GW170817 at merger time tm with accompanying GRB170817A at tGRB = tm + tg following a gap tg ≃ 1.7 s (bottom). The central engine of GRB170817A formed in this gap by causality, whereby tm < ts < tGRB for the start-time of accompanying gravitational radiation. This gap condition poses a discrete event time problem for the probability p of false alarm (PFA) for a candidate detection to be in the 1.7 s gap between GW1709817–GRB170817A: p = tg/T marked by the (unique) global maximum of an indicator function over an observational time T given the uniform prior on astrophysical event times (Appendix A). When sufficiently large, the energy output ℰGW marked by ts can break the degeneracy between neutron stars and black holes.

In the text

Fig. 2. Eight response curves of the two-level pipeline by χ-image analysis of merged (H1,L1)-spectrograms produced by butterfly matched filtering using a total of 56 signal injection experiments on LIGO O2 data covering GW170817 (top-left panel). These model injections are defined by Eq. (7). The dashed line is the mean of $\widehat{\chi }$-peaks averaged over (tiny) time-sides Δt within the light-travel time of 10 ms between H1 and L1 identified with extended emission to GW170817EE in a post-merger descending branch (lower panels). Precision realized is shown by parameter recovery (top-right panel). The parameter “field-of-view” covered for descending chirps covers energies 0 < ℰgw < 8% M⊙c2 (7 steps) and timescales of descent τs (eight steps, in units of seconds) with start frequency fs = 680 Hz for a model source distance D = 40 Mpc given by GW170817. These injection signals are in the shot-noise dominated regime of the LIGO detectors, where noise increases with frequency and hence sensitivity decreases with τs. Inserted is a plot of observed τs of post-merger emission to GW170817 (black square), overlaid to the calibration curve (blue curve) to infer uncertainty (along abscissa) from observed scatter (along the ordinate). Sample of a single injection of a combined merger and post-merger signal (1810–1825 s) in the snippet of H1L1-data containing GW170817 (tm = 1842.43 s) (reprinted from Paper II), shown in the lower panel. The present 56 injections Eq. (7) vary over energy ℰGW and time-scale of descent τs, keeping merger signal injection fixed at the parameters of GW170817.

In the text

Fig. 3. PDF of the start time ts (blue) of post-merger gravitational radiation to GW170817 following rejuvenation in gravitational collapse to a rotating black hole, setting the initial data to GRB170817A (top-left panel). The PDF shown (bin size 0.05 s) is superposed to estimates of the lifetime of the progenitor hypermassive neutron star formed in the immediate aftermath of the merger (reviewed in Murguia-Berthier et al. 2021). Highlighted is the 1.7 s gap (magenta) between merger time tm = 1842.43 s and the onset of GRB170817A. Post-merger emission to GW170817 appears as a descending chirp in the (H1,L1)-spectrogram merged by frequency coincidences |fH1 − fL1| < 10 Hz (lower panels). For the merged spectrogram, parameter estimation (Fig. 4) derives from time-slides less than the duration τ of the templates (11) against a background generated by larger time-slides. Post-merger gravitational radiation is identified with the formation of the central engine of GRB170817A, whose duration is constrained by the lifetime of black hole spin (bottom panel).

In the text

Fig. 4. Joint PDF(ts, τs) in Eq. (7) from merged (H1,L1)-spectrograms by χ-image analysis following time-slide analysis (top-left panels). Clustering is seen in the (ts, τs)-plane about a global maximum ${\widehat{\chi }}^{\ast }=8.76$ in the extended foreground over small time-slides Eq. (11). (The peak in the cluster in extended foreground (S1) is slightly larger than what is seen in conventional foreground (S0).) This cluster is extremal in the sense of Eq. (13). Tracing produces PDF’s of τs and the delay time ts − tm, where tm is the merger time of GW170817. There is notable consistency in τs with the duration $T_{90}^{8-70\mathrm{keV}}=2.9\pm 0.3$ s of GRB170817A (Pozanenko et al. 2018), expected from the timescale of black hole spin-down in interaction with surrounding high-density matter following gravitational collapse of the hyper-massive neutron star in the immediate aftermath of this merger. Clustering in τs and ts as a function of time-slide Δt in the post-merger emission to GW170817 extracted from (H1,L1)-spectrograms merged by frequency coincidences |fH1 − fL1| < 10 Hz (top-right panels). Contour plots of $\widehat{\chi }(t_s,\mathrm{\Delta }t)$ (maximal over τs, fs, f0) over all time-slides Δt ϵ [−32, 32] s (S0 − 2) in the bottom panels. The right panel zooms in on Δt ϵ [−1, 1] s, highlighting clustering within Δt ϵ [−τ, τ] set by the intermediate time scale τ ≃ 0.5 s of our template bank. Crucially, this cluster is extremal in the sense of Eq. (13) in extended foreground (S1).

In the text

Fig. 5. Sample of ($t_s^{*},{\tau }_s^{*}$) identified in χ-image analysis of individual H1- and L1-spectrograms (top-left panels), representing two elements from a total of nm = 320 H1- and L1-spectrograms (lower panels). Peak values indicated (white rectangles) have $\widehat{\chi }=7.0$ (H1) and $\widehat{\chi }=5.3$ (L1). Over the snippet of T = 2048 s of H1L1-data, time differences $\left[t_s\right]=t_{s,\mathrm{H}1}^{\ast }-t_{s,\mathrm{L}1}^{\ast }$ and $\left[{\tau }_s\right]={\tau }_{s,\mathrm{H}1}^{\ast }-{\tau }_{s,\mathrm{L}1}^{\ast }$ identified by local maxima of $\widehat{\chi }(t_s,{\tau }_s)$ over the 64 segments of 32 s and 320 trials of H1- and L1-spectrograms are, as expected, uncorrelated, evidenced by a correlation coefficient ρ[ts], [τs] = 1.29 × 10−4. (Right-top-panels) PDF1(ts) and PDF2(ts) of H1 and, respectively, L1 are centered around ts − tm = (0.92±0.09) s and τs = (3.00±0.09) s. Binning over δtres = 0.025 s, cross-correlation of the two PDFs (over each segment covering H1L1-data over T = 2048 s) shows a global maximum in segment 58 containing GW170817 consistent with zero time difference within δtres = 0.025 s.

In the text

Fig. 6. Chronicle (2) showing detailed event timing in Fig. 1. Highlighted is ℰGW by rejuvenation in gravitational collapse of the hyper-massive neutron star formed in the immediate aftermath of the merger at ≲ts in the 1.7 s gap G between GW170817–GRB170817A, in a Kerr black hole with spin-energy $E_J^{+}$ exceeding $E_J^{-}$ in its progenitor according to (3). $E_J^{+}$ amply accounts for ℰGW = (3.5 ± 1)% M⊙c2, breaking the degeneracy between neutron stars and black holes.

In the text

Fig. A.1. MC simulation on discrete event timing over n = 37 cells associated with a gap condition at cell 30, shown in magenta (top and middle panels). The event time, te, of a single run is the location of the maximum of an indicator function, χ(t), over all cells. A signal from a source in the gap (H1) increases the probability, p(te|H1), of detection at this location, above the PFA p(te|H0) = 1/n of stationary background noise (H0). The probability p(te|H1) increases monotonically with S/N, here derived from 1024 runs (lower panel).

In the text

Fig. B.1. (H1L1-detector spectral energy density Sn(f) for the T = 2048 s data snippet containing GW170817 following a band pass filter 10-1700 Hz (left panel). Violin modes appear prominently associated with the suspension of optics. Frequencies up to about 1700 Hz can be used in injection experiments. H1 and L1 detector noise is very similar though not identical during the GW170817 event. Whitening of H1 and L1 strain data by normalisation (B.1) produces flat spectra with lines suppressed (right panels).

In the text

Fig. C.1. Time-frequency domain illustrating filtering signals with frequencies slowly wandering in time, that appear as trajectories with finite time-derivative (top panel). Schematically, these trajectories passing through “butterflies" (green) defined by a finite slew rate (time rate-of-change) in frequency |df(t)/dt| ≥ δ > 0 (insert). It is realized by matched filtering over a dense bank of chirp-like time-symmetric templates of intermediate duration. (Reprinted from van Putten (2016).) Identification of broadband Kolmogorov spectrum in an ensemble averaged spectrum to the Nyquest frequency of 1 kHz of 42 bright long GRBs in the BeppoSAX catalog by butterfly filtering over a dense bank of 8.64 million templates, shown as a purple line (bottom panel). This result demonstrates sensitivity to turbulence better than an order of magnitude compared to conventional Fourier analysis (blue). Note: reprinted from van Putten et al. (2014a).

In the text

Fig. C.2. Master template serving as a seed to the construction of a template bank, given by an exponential decay in frequency f0 ≥ f(t) > f1 over 0 ≤ t ≤ T (top panel). Time-slicing produces chirps of intermediate duration P ≪ τ ≪ T, where P denotes the characteristic period of the (candidate) signal of interest. Time-symmetric templates obtain by linear superposition with their time-reverse (right), realizing equal sensitivity for ascending and descending chirp-like signals. (Reprinted from van Putten et al. (2014a).) Illustration of the suppression of constant frequency input (lower panel) in the form of the application of butterfly filtering to unwhitened LIGO data (left). Lines such as shown in Fig. B.1 fail to correlate with butterfly templates, leaving a smooth and approximately white background noise (right). Note: reprinted from van Putten (2016).

In the text

Fig. D.1. Performance overview of GPU-accelerated butterfly matched filtering on a quad-GPU node. Performance of clFFT in CSP/Out-of-Place as a function of array size N for two batch sizes measured by allocation M in global memory (left panel). Array sizes exceeding 212 cause a drop in performance as matrix transpose exceeds 32 kB of local memory, limited now by 1TB/s bandwidth to Global Memory in HBM2. We note that clButterfly uses a default size N = 217 (32 s of data at 4096Hz sampling rate) with batch size over n = 128 per frame of 4096 s. A quad-GPU node hereby features a maximum of about 5 teraFLOPs, circumventing limitations of PCIe bandwidth optimized by pre- and post-callback functions (van Putten 2017). The throughput of this quad-GPU node is shown up to six times real-time over a template bank normalized to size 220 (right panel). GPU-CPU communication passes through a post-callback function, limiting data-transfer to tails exceeding a threshold κ = 2 in normalized butterfly matched filtering output, effectively 0.1% of all matched filtering results. Throughput peaks in multi-frame analysis at a cross-correlation rate of about 200kHz.

In the text