The Central Engine of GRB170817A and the Energy Budget Issue: Kerr Black Hole versus Neutron Star in a Multi-Messenger Analysis

Abstract

Upcoming LIGO–Virgo–KAGRA (LVK) observational runs offer new opportunities to probe the central engines of extreme transient events. Cosmological gamma-ray bursts (GRBs) and core-collapse supernovae (CC-SNe), in particular, are believed to be powered by compact objects, i.e., a neutron star (NS) or black hole (BH). A principal distinction between an NS and BH is the energy reservoir in the angular momentum $E_J$. Per unit mass, this reaches a few percent in a rapidly rotating NS and tens of percent in a Kerr BH, respectively. Calorimetry by ${\mathcal{E}}_{GW}$ on a descending chirp may break the degeneracy between the two. We review this approach, anticipating new observational opportunities for planned LVK runs. GRB170817A is the first event revealing its central engine by a descending chirp in gravitational radiation. An accompanying energy output ${\mathcal{E}}_{GW}\simeq 3.5\%M_{\odot }{c}^2$ is observed during GRB170817A in the aftermath of the double neutron star merger GW170817. The progenitors of normal long GRBs, on the other hand, are the rare offspring of CC-SNe of type Ib/c. Yet, the extended emission to SGRBs (SGRBEEs) shares similar durations and the same Amati-relation of the prompt GRB emission of LGRBs, pointing to a common central engine. The central engine of these extreme transient events has, hitherto, eluded EM observations alone, even when including neutrino observations, as in SN1987A. The trigger signaling the birth of the compact object and the evolution powering these events is expected to be revealed by an accompanying GW signal, perhaps similar to that of GRB170817A. For GRB170817A, ${\mathcal{E}}_{GW}$ exceeds $E_J$ in the initial hyper-massive neutron star (HMNS) produced in the immediate aftermath of GW170817. It identifies the spin-down of a Kerr BH of mass ∼$2.4M_{\odot }$ defined by the total mass of GW170817. This observation is realized in spectrograms generated by Butterfly matched filtering, a time-symmetric analysis with equal sensitivity to ascending and descending chirps, calibrated by signal injection experiments. It is implemented on a heterogeneous computing platform with synaptic parallel processing in F90/C++/C99 under bash. A statistical significance of $5.5\sigma$ is derived from multi-messenger event timing, based on a probability of false alarm (PFA) factored over a probability $p_1=8.3\times {10}^{-4}$ by causality and a p-value $p_2=4.9\times {10}^{-5}$ of consistency between H1 and L1 observations. For upcoming observations, this approach may be applied to similar emissions from SNIb/c and GRBs in the Local Universe, upon the mass-scaling of present results by the mass of their putative black hole-central engines.

1. Introduction

With the detection of compact binary coalescence, which is well-established by LIGO–Virgo–KAGRA (LVK) [1,2,3], advancing search methods promise astronomical observations in a radically new window to unmodeled gravitational-wave (GW) transient sources. The evolving Universe is rife with extreme transient sources associated with the endpoints of massive stars and their compact stellar mass remnants, i.e., neutron stars (NS) and black holes (BH). Stellar mass BHs are believed to form in the core collapse of massive stars at the tail of the Salpeter stellar initial mass function [4]. Those observed with LIGO–Virgo–KAGRA (LVK) are in binaries sufficiently compact to merge within the Hubble time and appear to originate in the late-time Universe, no earlier than around the peak of the cosmic star formation rate [5]. LVK observations include a few mergers of double neutron star (DNS) and NS-BH systems. The DNS merger GW170817 stands out for its associated kilonova AT2017gfo and GRB170817A [6,7,8,9], providing the first member of the upcoming joint EM-GW surveys expected from upcoming LVK observational runs. With improved detector sensitivity, upcoming runs offer new opportunities to probe the formation and evolution of the central engines of core-collapse supernovae (CC-SNe) and gamma-ray bursts (GRBs).

Gamma-ray bursts were discovered serendipitously [10] as non-thermal high-energy transient events with a bimodal distribution in duration, isotropically distributed over the sky [11,12,13,14]. They appear as short ($T_{90}\lesssim 2$s) and long GRBs ($T_{90}\gtrsim 2$s), mirrored in relatively hard and soft extended spectra, respectively, with the latter satisfying the Amati-relation [15] between the peak energy ($E_p$) and isotropic equivalent energy output ($E_{iso}$), and related correlations [16,17,18,19,20,21]. While the prompt GRB emission has been attributed to dissipation in shocks derived from ultra-relativistic baryon-poor jets (BPJs) (e.g., [22,23,24]) and references therein), our description is mostly phenomenological and a detailed theoretical modeling is still lacking. Nevertheless, a BH appears favored over an NS in powering their central engines [25,26,27,28]. BPJs may be produced in a black hole–torus system [29,30,31,32] by virtue of baryon purity along open magnetic flux tubes subtended over the black hole event horizon [33,34].

GW170817 demonstrates a DNS origin of the short GRB GRB170817A [6,35], while normal long GRBs are known to be associated with the envelope-stripped CC-SNe of type Ib/c (e.g., [36]). The latter follows the identification of their cosmological origin [37,38,39,40,41] in the endpoint of massive stars [42,43,44,45,46,47,48]. Some short GRBs come with long-duration extended emissions (SGRBEEs), which satisfy the aforementioned Amati-relation of long GRBs [49,50]. GRB060614 stands out with a duration of $T_{90}=102$ s (15–350 keV) [51,52,53,54,55]. With no evidence of an associated supernova [51], its initial hard spike, characteristic of short GRBs, points to a common progenitor, believed to be a merger similar to GW170817, as opposed to a failed supernova. Attributing its initial short hard spike to the formation of its central engine, this would be obscured in a failed supernova by the opaque remnant stellar envelope. By similarity in duration and spectra of EE to SGRBs and the prompt GRB emission of LGRBs derived from mergers and, respectively, CC-SNe, these high-energy emissions appear to derive from a common central engine. Their long durations point to a secular time scale for the lifetime of their central engines, i.e., a dissipative rather than a dynamical time scale, such as the fall-back of matter from debris or, respectively, a remnant stellar envelope. This timescale possibly represents the lifetime of the spin of a rapidly spinning Kerr BH interacting with the surrounding high-density matter [56].

Corrected for beaming, the prompt GRB emission of LGRBs has true energy in gamma rays narrowly distributed about [57,58,59].

$$\begin{array}{c}\hfill E_{\gamma }\simeq 9\times {10}^{50}\mathrm{erg}.\end{array}$$

(1)

LGRBs show a positive correlation of luminosity with variability [60], the latter down to time scales $\delta t_{var}=\mathcal{O}$(ms) [25]. The extreme luminosity of GRBs and their possible counterparts in sources of ultra-high energy cosmic rays (UHECRs) may, in fact, originate from the intermittent behaviors of their central engines (e.g., [61,62,63,64,65]). The dimensionless measure of energetic time-variability [66,67]

$$\begin{array}{c}\hfill {\alpha }_E=\frac{E_{\gamma }}{L_0\delta t_{var}}=\mathcal{O}\left({10}^{-6}\right)\end{array}$$

(2)

far exceeds the same of Eddington-limited sources, such as the micro-quasar GRS1915+105 [68], showing their central engines to be extremely compact. Here, $L_0=c^5/G$ denotes the unit of luminosity in general relativity, given by the velocity of light c and Newton’s constant G. Combined, ms variability and (2) point to a central engine powered by an NS or stellar mass BH [33,69,70,71,72,73].

The energetic time-variability (2) ties in well with a stellar mass compact object produced in the core collapse of massive stars, when the nuclear burning phase ceases. Notably, pulsars in supernova remnants are neutron stars, e.g., Vela and the Crab [74], and some may be found in NS binaries, such as the Hulse–Taylor system PSR 1913+16 [75] and the relativistic binary PSR J0737-3039 [76].

While the identification of pulsars with rotating neutron stars is unambiguous, the electromagnetic spectrum offers compelling but inconclusive evidence of the nature of the central engine of GRBs despite their common origin in (the remnants of) massive stars. From the first principles, a BH cannot act similar to a pulsar by the inherent alignment of the magnetic moment and spin [77,78]. On account of the Kerr metric [79], however, the energy reservoir per unit mass in the angular momentum $E_J/M$ of a rotating BH can vastly “outsize” that of a rapidly spinning NS by an order of magnitude. This offers the potential to break the degeneracy between the two, provided a suitable fraction of $E_J$ is channeled into radiation “unseen”—beyond the electromagnetic spectrum—given the existing constraint (1).

LVK observations offer, for the first time, a means to extend our probe of central engines to all current radiation channels, including gravitational radiation.

Searches for unmodeled GW emissions may be approached in spectrograms, covering the bandwidth of sensitivity of the LVK detectors. Candidate GW transient emissions may, herein, appear as chirps, ascending from mergers or descending in the spin-down of a compact object, notably a neutron star or a black hole.

This appears opportune for extreme transient events powered by a BH interacting with surrounding high-density matter, when gravitational radiation may be the dominant channel of emission with the accompanying output in MeV neutrinos and electromagnetic radiation [34]. Starting from the progenitor, we look for the trigger signaling the birth of the central engine. Of particular interest is the DNS merger GW170817 [6] and, perhaps, the more numerous energetic SNIb/c associated with LGRBs [80,81].

At the relative proximity of $D\simeq 40$ Mpc [82,83], GW170817 produces a compact object triggering the onset of GRB170817A [6,7,8,9]. Extended emissions ${\mathcal{E}}_{GW}$ in gravitational radiation may tell us about the nature of the central engine of GRB170817A—NS or BH—in addition to the physical properties of the progenitor GW170817 [84,85,86,87], realizing a key objective of multi-messenger observations [88,89].

GW170817 reveals the extended evolution in complex multi-messenger emission, which includes the kilonova AT2017gfo and GRB170817A of duration $T_{90}^{8-70\mathrm{keV}}=\left(2.9\pm 0.3\right)$s [90], following this DNS merger at time $t_m$ [91],

$$\begin{array}{c}\hfill \begin{array}{ccccc}\mathrm{NS}+\mathrm{NS}\hfill & \to \hfill & \mathrm{HMNS}\hfill & \to \hfill & \mathrm{BH}+\mathrm{AT}2017\mathrm{gfo}+\mathrm{GRB}170817\mathrm{A}+\cdots ,\hfill \end{array}\end{array}$$

(3)

where dots refer to emissions in MeV neutrinos and gravitational radiation. AT2017gfo has been attributed to a powerful MeV neutrino-driven wind from the initial hypermassive NS (HMNS) produced in the immediate aftermath of GW170817, although additional contributions following delayed gravitational collapse cannot be excluded by (neutrino-driven) neutron-rich winds from a relatively lower temperature torus [33,34] around a BH [92,93,94,95,96].

A LIGO gravitational-wave observation of (3) requires a sensitivity close to the limit defined by the strain-noise amplitude of the detectors conforming to the theory of matched filtering [97]. In particular, the ascending chirp GW170817 may be followed by an unmodeled descending chirp characteristic for the spin-down of the HMNS [98,99,100] and/or the compact object powering GRB170817A. A suitable starting point to probe for unmodeled signals in (3) is provided by the Butterfly matched filtering over a bank of time-symmetric templates of intermediate duration, originally developed to identify a broadband Kolmogorov spectrum in bright GRBs of the BeppoSAX catalog at high frequencies [101]. This approach realizes equal sensitivity to ascending and descending chirps at a sensitivity on par with CBC when this bank is sufficiently dense in $(f,\dot{f})$: the frequency f and its time rate-of-change $\dot{f}$, suitably accelerated by heterogeneous computing [102].

We quantify the above-mentioned trigger by event timing. Given the uniform prior on timing in astrophysics on time scales less than a Hubble time and in instrumentation, event timing produces estimates of significance expressed in terms of the probability of false alarm (PFA): a type I error, given the null hypothesis of the background, assuming no event. Accurate event-timing parameterization is derived from probability density functions (PDFs), particularly for the start times $t_s$ of burst events. Combining statistically independent PFAs requires distinguishing between continuous and discrete stochastic variables (e.g., [70,103,104,105,106]). This is opportune when event timing extends over the mean and the difference of observations across two detectors, notably LIGO H1 and L1. Equivalent to a unitary transformation given by the rotation over $\pi /4$, this generates statistically independent PFAs. A leading-order interpretation of results derives from energy ${\mathcal{E}}_{GW}$ and frequency in the gravitational-wave emission at hand, subject to priors, such as the total mass-energy of the system and limits on the mathematical–physical properties of HMNS and Kerr BHs.

In this review, we revisit a recent direct detection of a descending chirp in gravitational waves during GRB170817A and the observational consequences for its central engine. Given the upcoming LVK observations, we emphasize the methods of analysis by anticipating new opportunities to probe the central engines of GRBs and associated SNIb/c in the Local Universe, given the significant improvement in sensitivity in upcoming LVK runs. Accordingly, this review is organized as follows:

• The spin energy in the gravitational collapse of the initial HMNS (e.g., [100,107,108]), to a Kerr black hole (Section 2);
• Event timing in multi-messenger observations (Section 3);
• Observing GW transient sources in time-symmetric spectrograms (Section 4);
• What triggered GRB170817A following GW170817? (Section 5);
• High-resolution event timing produced by exascale computing (Section 6, Section 7 and Section 8).

Section 7 and Section 8 are included for completeness, closely following the original work in [109].

In Section 9, we conclude with a discussion on observational consequences and an outlook on possibly similar approaches to study the central engine of GRBs and SNIb/c in the Local Universe with upcoming LVK observations in joint EM-GW observations.

2. Gravitational Collapse of an HMNS to a Kerr Black Hole

The multi-messenger events SN1987A and GRB170817A both have known progenitors: the blue supergiant (Sanduleak-69 202) in the LMC at a distance of 50 kpc [110,111] with an associated MeV neutrino burst [112,113] and, respectively, from the merger event GW170817 of a double neutron star (DNS). They have the initial formation of a (hyper)massive NS in common, whose fate is inconclusive based on these progenitor data alone. For a complete interpretation, it is essential, however, to determine its continuing evolution. (This includes its potential prior evolution in the case of SN1987A, below). A continued (delayed) gravitational collapse to a BH may have dramatic implications for the central engine powering the supernova explosion and, respectively, GRB170817A, given the principal distinction in the energy reservoir in the angular momentum $E_J$ of an HMNS [114] and a Kerr BH (e.g., [115]). After all, classical mechanics teaches us that the prompt gravitational collapse increases $E_J$, preserving mass M and angular momentum J.

In (3), $E_J$ gradually decreases during the merger phase GW170817, reaching the first minimum at coalescence to an HMNS. If long-lived, $E_J$ may continue a decline upon radiating baryon-rich winds and, perhaps, gravitational radiation when sufficiently non-axisymmetric deformations are present. Any loss of stability, however, by spin-down, cooling, and internal re-distribution of J, may result in the delayed gravitational collapse to a BH. If so, this enhances $E_J$ as the following estimate shows. Conserving M and J in the prompt gravitational collapse preserves the specific angular momentum $a=J/M$ of the newly formed black hole and, hence, the dimensionless Kerr parameter $a/M=sin\lambda$. The spin energy $E_J^{-}\simeq (1/2)J\mathsf{\Omega }$ of the HMNS increases to $E_J^{+}=2M{c}^2{sin}^2(\lambda /4)$ [56,79]. Specialized to GW170817, we estimate [109]

$$\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}$$

(4)

where $J=I\mathsf{\Omega }$ for the moment of inertia $I\simeq (2/5)M{R}^2$ of the HMNS with the gravitational radius $R_g=GM/c^2$, radius R, and angular velocity $\mathsf{\Omega }$. Here, G is Newton’s constant and c is the velocity of light, and we ignore further enhancement from, e.g., hyperaccretion (e.g., [116,117]). This enhancement sets the stage for potentially significant extended emissions ${\mathcal{E}}_{GW}$ by coupling to non-axisymmetric high-density matter surrounding the black hole at relatively moderate gravitational-wave frequencies [34], distinct from what might derive directly from the progenitor HMNS.

Previously, we reported on the direct detection of a descending chirp in gravitational radiation during GRB170817A [118]. Starting between the time of merger and the onset of GRB170817A, its total energy output [81]

$$\begin{array}{c}\hfill {\mathcal{E}}_{gw}=(3.5\pm 1)\%M_{\odot }{c}^2\simeq 6.3\times {10}^{52}\mathrm{erg}\end{array}$$

(5)

is emitted at gravitational-wave frequencies $f_{gw}\lesssim 700$Hz. For quadrupole gravitational-wave emissions, (5) is to be seen against the estimated bound on $E_J^{-}$ in the spin energy of the initial HMNS,

$$\begin{array}{c}\hfill E_J^{-}=\frac{{\pi }^2}{5}{f}_{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}$$

(6)

With (5) exceeding (6) by a factor of about 4, even at a maximal radius R, Equation (5) breaks the degeneracy between the HMNS and a BH. This holds true especially given the start time $t_s\lesssim 1$s of (5), the post-merger to GW170817, when neutrino cooling and internal dissipation are expected to produce a uniform rotator [119]. In contrast, Equation (5) is consistent with $E_J$ radiated off in the spin-down of a Kerr BH against a non-axisymmetric high-density torus [81].

The above shows the power of calorimetry in gravitational radiation, breaking the degeneracy between an (HM)NS and BH formed in continuing gravitational collapse. Indeed, EM calorimetry falls woefully short; the combined output of GRB170817A and AT2017gfo [120,121,122,123,124] is limited to

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

(7)

Even so, EM observations provide us with crucial information on event timing.

3. Event Timing in Multi-Messenger Observations

SN1987A is instructive once again, for its similar arrival times of the MeV neutrino burst in Kamiokande II and IMB [112,125,126], signaling the core collapse of its massive stellar progenitor, well before the optical light of this type II supernova explosion in the LMC at a distance of about 50 kpc. Likewise, GW170817 is picked up by LIGO H1 and L1 detectors, essentially simultaneously (L1 slightly before H1) to the progenitor of GRB170817A in NGC4993 at a distance of 40 Mpc. (Accompanying the MeV neutrino emission is a posteriori undetectable). These two delay times represent causality between the progenitor star, trigger, and the ensuing events: SN 1987A and, respectively, the binary progenitor of the ensuing high-energy transient GRB170817A.

While consistent, the two multi-messenger datasets in SN1987A and GRB170817A do not necessarily cover the complete evolution from the progenitor to the final remnant. Moreover, Equation (3) shows an additional kilonova AT2017gfo, derived from neutron-rich winds in evolution after the merger, which may include the delayed gravitational collapse to and evolution of a rotating black hole. SN1987A might have had two precursor triggers in the MeV emission (see [113]). A low-count but potentially significant MeV neutrino burst was picked up by the Liquid Scintillation Detector (LSD) and KII a few hours prior to the main KII-IMB MeV neutrino burst (reviewed in [127,128]). It may suggest a possible two-stage process, starting with the formation of a DNS to radiate off the excess angular momentum in the gravitational radiation before coalescence into a (proto)-NS producing the primary MeV neutrino burst [126].

3.1. The Gap Time $t_g$ between GRB170817A and Its Progenitor GW170817

A complete observation covering the entire evolution from the progenitor to the final remnant in (3), therefore, may involve additional precursor GW signals. In particular, the true trigger of GRB170817A may be defined by a gravitational-wave signal that starts in the gap between the merger GW170817 at $t_m$ and the time-of-onset of GRB170817A at $t_{GRB}$,

$$\begin{array}{c}\hfill t_g=t_{GRB}-t_m\simeq 1.7\mathrm{s}.\end{array}$$

(8)

This gap time may be accounted for by the sum of the lifetime of the HMNS on the basis of the duration $t_w$ of the MeV neutrino wind powering AT2017gfo [129] and a delay time $\Delta t_{GRB}$ between the observed time-of-onset $t_{GRB}$ of GRB170817A,

$$\begin{array}{c}\hfill t_g=t_w+\Delta t_{GRB}.\end{array}$$

(9)

Following the gravitational collapse of the HMNS, we identify $\Delta t_{GRB}$ with the relativistic delay time between GRB170817A and the start time $t_s\simeq t_m+t_w$ of the gravitational-wave emission signaling the birth of its engine. Attributing GRB170817A to a dissipative front at a distance r from the central engine, propagating at an angle $\theta$ relative to the line-of-sight to the observer, $\Delta t_{GRB}=\left(1-\beta cos\theta \right)r/c$; thus,

$$\begin{array}{c}\hfill r=\frac{c\Delta t_{GRB}}{\left(1-\beta cos\theta \right)}\le c\Delta t_{GRB}\times min\left[2{\Gamma }^2,\frac{1}{1-cos\theta }\right],\end{array}$$

(10)

assuming the same velocity of propagation for both gravitational waves and electromagnetic radiation (e.g., [130,131]) (also see [22]).

Assuming $\theta \lesssim 1/\Gamma$ for a Lorentz factor $\Gamma$ of the ultra-relativistic jet powering the GRB, $\Delta t_{GRB}\simeq r/2c{\Gamma }^2$ (cf. [23,132]). However, GRB170817A is likely seen off-axis at $\theta \simeq {20}^o$[133]. Estimating r, this approximation may not apply, and the more conservative second term on the right of (10) must be used instead.

Assuming AT2017gfo to be mostly powered by the MeV neutrino-driven neutron-rich wind from the initial HMNS, its “blue” component [134,135] implies a finite lifetime of this HMNS of about $t_w\lesssim 1$s [129,136,137,138]. (see [139,140,141,142] for a possibly long-lived supermassive NS remnant). This still leaves a positive delay time $0<t_g-t_w\lesssim 1$ between the gravitational collapse of the HMNS and GRB170817A, here identified with the relativistic time delay (10). We proceed with the time of formation of the central engine of GRB170817A, i.e., at a preceding trigger time $t_s$,

$$\begin{array}{c}\hfill C_1:t_s-t_mϵG\end{array}$$

(11)

where $G=\left[t_m,t_{GRB}\right]$ denotes the interval of (8) (Figure 1) with an observable descending chirp (Figure 2, Section 5) starting at

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

(12)

based on the mean of $t_s$, inferred from an independent analysis of H1 and L1 data (Section 8). This delay time of the gravitational collapse (12) is indeed consistent with $t_w=\left(0.98\pm 0.3\right)$ s [129] of the HMNS inferred from AT2017gfo [24,81,129,143,144,145,146,147,148] (reviewed in [9]).

Combined with $\Gamma \simeq 70$ [133], $\Delta t_{GRB}=t_g-t_s\simeq 0.8$ s, we arrive at (cf. [149])

$$\begin{array}{c}\hfill 0.03\mathrm{AU}\lesssim r\ll 16\mathrm{AU}\end{array}$$

(13)

for $1\ll \theta \lesssim 20°$, where AU $=1.5\times {10}^{13}$ cm. For the first time, EM-GW observations hereby constrain (10), which is central to the process of prompt GRB emissions.

3.2. Event Timing by the Trigger Time $t_s$ of GRB170817A

In searches for gravitational-wave emissions triggered by gravitational collapse, Equation (11) defines the aforementioned gap condition, providing a causality constraint on $t_s$ with a uniform prior. These priors are considered canonical for timing in astrophysics, (specifically over time scales less than a Hubble time) and are also widely adopted in instrumentation. Over an observational time T, Equation (11) hereby carries PFA${}_1$—a type I error—for it to be satisfied by mere chance under the null hypothesis ($H0$) of the stationary detector noise,

$$\begin{array}{c}\hfill p_1=\frac{t_g}{T}.\end{array}$$

(14)

A specific implementation involves an indicator function $\chi \left(t\right)$ that marks an event time $t_e$, here, the start time $t_s$ of a candidate’s emission feature, according to a global maximum of $\chi \left(t\right)$ over the observation time T. $\chi \left(t_e\right)$ is a continuous stochastic observable, where $t_e$ is well-defined and unique. Under the null hypothesis, $t_e$ is uniformly distributed over $[0, T]$, and (11) is satisfied by mere chance with probability $p\left(t_e\right|H0)=1/n$, following a partitioning of $[0, T]$ over $n=T/t_g$ bins. This produces PFA (14).

In contrast, gravitational-wave emissions from a putative central engine of GRB170817A ($H1$) change our expectation to $p\left(t_e\right|H1)$, with a bias to (11), provided that $\chi \left(t\right)$ is suitably devised to pick up its emission. Accordingly,

$$\begin{array}{c}\hfill p\left(t_e\right|H1)\ge p\left(t_e\right|H0)\end{array}$$

(15)

for $t_eϵG$, depending on the signal-to-noise ratio (SNR). A Monte Carlo simulation of Pr$\left(t_eϵG\right)$ illustrates a function of the SNR, where $t_e$ is the location of the maximum of $\chi \left(t\right)$ in data over n bins given by the sum of noise plus a signal from a source at G in a selected bin [109]. In the absence of a signal, this probability reduces to the PFA, which, here, is equal to the chance of winning a single bet in European roulette (e.g., [150]).

Event timing by (11) produces a PFA (14) fixed by context; here, it is GW170817-GRB170817A. This is quite different from PFAs inferred from a statistic in a continuous observable, such as SNR. In particular, (14) does not involve numerical or observational parameters, such as trial factors, provided that (11) applies. The latter requires that the measurement uncertainty in $t_e$ be small relative to $t_g$, i.e.,

$$\begin{array}{c}\hfill {\sigma }_{t_s-t_m}\ll t_g\ll T,\end{array}$$

(16)

where ${\sigma }_{t_s-t_m}\simeq {\sigma }_{t_s}$ is derived from the PDF($t_s$) defined by the measurement at hand (Section 7). In this sense, (11) represents a discrete event-timing condition. The minimum of Pr$(t_s-t_mϵG)$, PFA $p_1$ (14) is an ordinary probability that can be easily combined with other PFAs [151,152,153]. Box 1 highlights some key elements of this section.

Box 1. The delayed trigger time $t_s$ of GRB170817A.

• GW170817 initially produced an HMNS, evidenced by the kilonova AT2017gfo, followed by GRB170817A across a duration gap of $t_g=t_{GRB}-t_m\simeq 1.7$s, containing the time of the gravitational collapse of the HMNS to a BH, $t_s$.
• The black hole is initially rapidly spinning, as evidenced by an output ${\mathcal{E}}_{GW}\simeq 3.5\%M_{\odot }{c}^2$ in $f_{GW}\lesssim 700$Hz, exceeding the spin energy $E_J^{-}$ of the HMNS.
• The EM-GW event timing shows tight consistency in the lifetime $t_w=\left(0.98\pm 0.3\right)$s of the HMNS [129] and the start time $t_s-t_m=\left(0.92\pm 0.083\right)$s of the extended emission (Figure 2).
• Based on EM-GW data, a relativistic time delay (10) shows that GRB170817A is derived from a dissipative front at a distance (13) from the central engine.

4. Observations of GW Transient Emission in Time-Symmetric Spectrograms

The prototypical spectrograms are produced by time-sliced fast Fourier transforms (FFTs). High-resolution spectrograms are obtained by applying a filtering process that covers the effectively dense cover of the desired bandwidth by discrete frequencies, subject to the Nyquist sample rate criterion. When applied to searches for transient signals in noisy data [97,154], a further restriction is derived from the choice of the time slice $\tau$: a finite bandwidth $\Delta f\simeq 1/\tau$ associated with each frequency is derived from the time–frequency uncertainty relation,

$$\begin{array}{c}\hfill \tau \Delta f\gtrsim 1.\end{array}$$

(17)

This sensitivity allows for the detection of transient signals characterized by time-varying frequencies $f\left(t\right)$ with a finite slope $\dot{f}=df\left(t\right)/dt$ up to $\Delta f/\tau \simeq 1/{\tau }^2$. Time-sliced FFT hereby realizes a time-symmetric sensitivity, improved by a factor $\sqrt{N}$ with $N\simeq f\tau$, according to the theory of matched filtering. When $\dot{f}$ is large, the condition $\tau \lesssim 1/\sqrt{\left|\dot{f}\right|}$ limits $\tau$. According to the theory of matched filtering, the signal-to-noise improvement is hereby limited to

$$\begin{array}{c}\hfill \sqrt{N}\lesssim f^{\frac{1}{2}}{\left|\dot{f}\right|}^{-1/4}.\end{array}$$

(18)

Figure 3 shows the spectrograms by the time-sliced FFT to H1- and L1-data of GW170817.

For instance, a chirp-changing frequency by a few hundred Hz over the course of a few seconds has a slew rate of $\dot{f}\simeq 100$s${}^{-2}$ and, hence, $\tau \lesssim 0.1$. For frequencies below 1 kHz, (18) leaves an improvement of no more than a factor of a few.

4.1. Butterfly Matched Filtering

To go beyond (18), we circumvent (17) by considering an alternative filter bank of chirp-like templates of frequencies with finite $\dot{f}$, with compact support over a choice of $\tau$. Chirp-like templates enable the capture of signals with finite $\dot{f}$ over some intermediate time scale $\tau$ of the phase coherence. (In searches for unknown signals, $\tau$ is a search parameter). To this end, we consider a two-dimensional bank of templates densely covering a region in $(f,\dot{f})$. Time symmetry can be preserved by construction, i.e., by taking the differences between such chirp-like templates and their time-reversed counterparts. This approach has been successfully applied in the identification of the broadband Kolmogorov turbulence in the Poisson noise-dominated light curves of bright GRBs in the BeppoSAX catalog [101]. Figure 4 shows the resulting spectrum, which can be resolved with an order of magnitude beyond what is achieved with the conventional Fourier analysis. An early porting to the gravitational-wave data analysis to SN2010br is reported in [155].

In searches for unknown signals, a key parameter is the intermediate time scale of the phase coherence $\tau$, which may be present in the signal with a minimum $\delta >0$ in the time ratio-of-change of frequency (Figure 4). The core idea of Butterfly matched filtering involves using a filter composed of a bank of time-symmetric templates with a duration $\tau$; these templates are dense in both f and $\dot{f}$. $\tau$ should be restricted to the values intermediate to the total duration of the candidate signal and the total search period. Additionally, the change of frequency in time satisfies

$$\begin{array}{c}\hfill \left|\frac{df\left(t\right)}{dt}\right|\ge \delta >0.\end{array}$$

(19)

Figure 4 provides a schematic representation using butterflies, indicating a finite range of slew rates (19) at each point in the time–frequency domain.

Time-symmetric templates are created by linearly superposing a segment from a master template with its time-reversed counterpart. The resulting time-symmetric templates ensure equal sensitivity to ascending and descending chirps [101],

$$\begin{array}{c}\hfill u_k\left(t^{\prime }\right)=z_k\left(t^{\prime }\right)-z_k(\tau -t^{\prime })(0\le t^{\prime }\le \tau ),\end{array}$$

(20)

where $z_k\left(t^{\prime }\right)$ refers to a chirp-like template with compact support over a duration $\tau$.

Generally, the performance gain in picking up signals with time-varying frequencies via matched filtering is expected to depend on the type of noise. The Poisson noise-dominated analysis of Figure 4 shows an improvement in sensitivity by (20), up to an order of magnitude over the conventional Fourier analysis. The gravitational strain-noise data can be effectively modeled as colored Gaussian noise, and further, with numerous lines, notably the violin modes of suspension wires. While Poisson noise is tightly correlated with the signal strength, the gravitational-wave detector noise is completely decoupled from signals passing by. The sensitivity gain of Butterfly matched filtering over the Fourier analysis may, therefore, be similar but not identical to the result shown in Figure 4. An explicit calculation is obtained by appealing to the theoretical gain in the signal-to-noise ratio by $\sqrt{N}$ in the matched filtering, where N denotes the number of periods captured phase-coherently over the intermediate time scale $\tau$ set by the duration of a template.

4.1.1. Sensitivity Gain over the Time-Sliced Fourier Analysis

Regarding signals with time-varying frequencies, the gain in sensitivity of Butterfly matched filtering over the time-sliced Fourier analysis can be seen as follows. Consider a quasi-monochromatic signal admitting a Taylor series expansion $f\left(t\right)=f_0+t{f}_1+\frac{1}{2}{t}^2{f}_2+\mathcal{O}\left(t^3\right)$, where $f_n=d^{n}f\left(t\right)/d{t}^n$ refers to the n-th derivative at $t=0$ of its frequency $f\left(t\right)$ as a function of time t. The matched filtering over a duration $\tau$ attempts to track the associated change in phase, $\phi \left(\tau \right)=2\pi \left(f_0\tau +\frac{1}{2}{f}_1{\tau }^2+\frac{1}{6}{f}_2{\tau }^3+\mathcal{O}\left({\tau }^4\right)\right)$. The Fourier analysis captures ${\phi }_0=2\pi f_0\tau$ over a duration $\tau$ until this is out of phase by $\pi$, relative to the true total phase $\phi \left(\tau \right)$, causing destructive interference. That is, $\frac{1}{2}{f}_1{\tau }^2\lesssim 1/2$, i.e., $\tau \lesssim {\tau }_A$, ${\tau }_A=1/\sqrt{f_1}$ is commensurate with the above-mentioned derivation based on the time–frequency uncertainty relation. Extending this argument to Butterfly matched filtering, which captures ${\phi }_1=2\pi \left(f_0\tau +\frac{1}{2}{f}_1{\tau }^2\right)$, we insist that $\frac{1}{6}{f}_2{\tau }^3\lesssim 1/2$, i.e., $\tau \lesssim {\tau }_B$ with ${\tau }_B={\left(3/f_2\right)}^{1/3}$. The improvement in $N\simeq f_0\tau$ is given by the ratio of these bounds, ${\tau }_B/{\tau }_A=3^{1/3}{f}_1^{1/2}/f_2^{1/3}$, implying a sensitivity improved by

$$\begin{array}{c}\hfill \mathrm{gain}\simeq 3^{1/6}{f}_1^{1/4}{f}_2^{-1/6}\end{array}$$

(21)

of Butterfly matched filtering over the time-sliced Fourier analysis.

4.1.2. Application to Gravitational-Wave Data

By $\delta >0$ in (19) being sufficiently large, constant frequency signals hereby tend to be relatively suppressed, distinguishing Butterfly matched filtering from conventional Fourier analysis. This effective suppression of constant frequency is demonstrated in the line suppression of unwhitened LIGO data [155], and is further improved by pre-whitening [109,118].

An illustrative example is a signal featuring a change in frequency $\Delta f$ with the characteristic time scale of decay ${\tau }_s$. In this case, $f_1\simeq \Delta f/{\tau }_s$ and $f_2\simeq \Delta f/{\tau }_s^2$ and, hence,

$$\begin{array}{c}\hfill \mathrm{gain}\simeq 3^{1/6}{\left(\Delta N\right)}^{1/12},\end{array}$$

(22)

where $\Delta N={\tau }_s\Delta f$. By applying (22) to a descending chirp in gravitational radiation below 700 kHz down to 200 Hz over ${\tau }_s=3$s, we have $\Delta f=500$Hz and $\Delta N=1500$. This leads to a gain of $\simeq$2.2. In the application to gravitational-wave observations, this gain effectively brings the source closer by a factor of 2. This factor of 2 represents a significant improvement in search sensitivity roughly equivalent to one decade improvement in LIGO sensitivity.

The results of the individual H1 and L1 spectrograms may further be considered in (H1,L1)-spectrograms merged by frequency coincidences (Figure 5), providing further improvement, in contrast, by a factor of $\sqrt{2}$. The total gain over the time-sliced Fourier analysis of a single detector output is, hereby, 3.12.

By extending (21) and (22) to include dense filtering up to $f_n$, we expect a gain over the time-sliced FFT:

$$\mathrm{gain}=\sqrt{{\displaystyle \frac{{\tau }_n}{{\tau }_0}}}\simeq f_1^{\frac{1}{4}}{\left({\displaystyle \frac{(n+2)!}{2f_{(n+1)}}}\right)}^{\frac{1}{2(n+1)}}\simeq {\left({\displaystyle \frac{(n+2)!}{2}}\Delta N^{\frac{n}{2}}\right)}^{\frac{1}{2(n+2)}}\simeq \left\{\begin{array}{c}1n=0\hfill \\ 2.2n=1\hfill \\ 3.4n=2\hfill \end{array}\right..$$

(23)

For instance, when using $n=2$ and dense filtering in $(f,\dot{f},\ddot{f})$, an additional gain of $\simeq$1.54 over Butterfly matched filtering ($n=1$) is achieved, as well as a gain of $\simeq$4.80 over time-sliced FFT ($n=0$) in the merged (H1,L1) spectrograms. Indeed, Figure 5 shows improved visual contrast consistent with the gain (23) for the ascending–descending chirp in Figure 3. Box 2 highlights some key findings of this section.

Box 2. High-resolution spectrograms: Fourier versus Butterfly matched filtering.

• GW transients may be ascending or descending chirps in gravitational radiation of mergers and, respectively, the spin-down of a compact object.
• GW transients can be searched for in high-resolution spectrograms covering the bandwidth of sensitivity of LVK over a time scale $\tau$ of time-slicing.
• In Fourier-based spectrograms, sensitivity (18) is limited by the time–frequency relation (17).
• A gain in sensitivity (23) obtains in Butterfly matched filtering over a bank of time-symmetric chirp-like templates densely covering a region of $\left(f,\dot{f}\right)$ spaces.

5. Ascending–Descending GW Transient Emission during GW170817-GRB170817A

Quite generally, CBC and remnant compact objects shedding angular momentum in gravitational radiation are manifested in ascending and, respectively, descending chirps. Butterfly matched filtering [155] aims to capture ascending and descending chirps alike, selectively amplifying signals whose frequencies are slowly wandering in time and suppressing signals whose frequencies are essentially constant.

5.1. Observation in H1L1 Data

The result (5) is seen in (H1,L1)-spectrograms generated by Butterfly matched filtering. It has equal sensitivity to ascending and descending chirps by construction (20), all else being equal. Previous LIGO searches by power-excess methods, as discussed in [6,156,157], have been found to fall short of a detection threshold of ${\mathcal{E}}_{gw}\gtrsim 650\%M_{\odot }{c}^2$ [158]. This search threshold is within the excluded zone [109], beyond the total mass-energy of GW170817, with a sensitivity to signals post-merger being 0.3% compared to ${\mathcal{E}}_{gw}=2.25\%M_{\odot }{c}^2$ in GW170817 [6].

Figure 5 shows a probe of H1L1 data covering GW170817-GRB170817A by Butterfly matched filtering. The extended ascending–descending emission suggests gravitational radiation contemporaneous with GRB170817A. Produced by its central engine, this would represent emissions from compact remnant(s) of GW170817.

For a detailed consideration of the extended emission feature in Figure 5, we calibrate the sensitivity of our spectrograms to estimate the energy output ${\mathcal{E}}_{GW}$ in the descending branch as a potential constraint on the nature of the central engine of GRB170817A, followed by the parameter estimation and estimates of statistical significance.

5.2. Calibrated Response Curves

Sensitivity can be calibrated using response curves generated through signal injection experiments. Here, we focus on signal characteristics for the spin-down of a compact object: a descending chirp in gravitational radiation (5).

For descending chirps resulting from the spin-down of a compact object due to gravitational radiation, the following “1 + 3” parameterization is used:

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

(24)

appended as a delayed descending chirp to a model ascending chirp of GW170817. These curves comprise a start time $t_s$—the trigger time—and three nuisance parameters, the $({\tau }_s,f_s,f_0)$ time scale of descent ${\tau }_s$, the start frequency $f_s$, and the late-time frequency $f_0$. Alternatively, (24) is “2 + 2” when estimating $(t_s,{\tau }_s)$ with nuisance parameters $(f_s,f_0)$.

Figure 6 shows the results of the descending chirp models following a model merger signal GW170817 at its distance $D=40$Mpc. These injections are performed before GW170817A, ensuring a sufficient distance to avoid interference with GW170817, yet remaining close enough to essentially share the same data quality.

Thus, Figure 6 shows two ascending–descending signals: the first corresponds to the injected signals, and the second represents GW170817 along with its accompanying delayed extended emission feature. Figure 6 includes results of the two-parameter injection experiments over energies and time scales

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

(25)

given the source distance $D=40$ Mpc of GW170817 [83], keeping $f_s=680$ Hz, and $f_0=98$ Hz.

A descending chirp of the form (24) is captured by scanning over all the parameters $(t_s,{\tau }_s,f_s,f_0)$, by $\chi$ image analysis over the H1L1 data of $T=2048$ s, containing GW170817. This process is realized in $2^{14}$ steps in (${\tau }_s,f_s,f_0$) at 33 Hz in $t_s$ for a total of ${10}^9$ values [109]. The response curves of Figure 6 are averaged over 41 tiny (surrogate) time slides between H1 and L1, which corresponds to their light crossing travel time of 10 ms. Box 3 some key conclusions on the central engine of GRB170817A.

Box 3. What is the central engine of GRB170817A?

• The ascending–descending chirp in GW170817 represents a merger followed by delayed spin-down of a Kerr black hole during GRB170817A.
• This observation is seen in spectrograms generated by Butterfly matched filtering calibrated by signal injection experiments.
• EM-GW event timing shows consistency between $T_{90}^{8-70\mathrm{keV}}=\left(2.9\pm 0.3\right)$s of GRB170817A [90], $T_{GW}\simeq 3.7$s and time scale of descent ${\tau }_s=\left(3.0\pm 0.1\right)$s.

6. Exascale Computing by Synaptic Parallel Processing

Observations of unknown signals require two (or more) detectors with comparable nominal behaviors to cross-correlate independent candidate signals. Over large datasets, such as full observational runs, the generation of spectrograms through matched filtering presents a novel computational challenge in generating spectrograms produced by matched filtering, covering a sufficiently broad range and density in $(f,\dot{f})$. For a candidate feature in a spectrogram, the next step is to quantify its properties by parameter estimation. This notably includes the extraction of the probability density function (PDF) of a choice of parameters, such as those encountered in event timing. Matched filtering is a computer-intense task involving cross-correlation between the data and templates, which scales bilinearly with the data length and size of the template bank. Substantial additional computational costs are, therefore, incurred, particularly in the process of parameter estimation.

6.1. Acceleration by High-Performance Computing

To meet this challenge, we appeal to modern computational approaches in hardware and software. Specific methods to accelerate matched filtering by a few orders of magnitude derive from efficient algorithms, modern computing hardware, and network computing:

• Efficient evaluation in the Fourier domain using the fast Fourier transform (FFT);
• Heterogeneous computing by offloading inverse FFTs to graphics processor units (GPUs);
• Distributed computing on a platform, load-balanced by synaptic parallel processing.

The requirement for these computational components is more evident by illustrating the scale of computation in obtaining a single spectrogram, using cross-correlation, for instance, for H1L1 data around GW170817-GRB170817A. Consider data with durations of $T=2048$ s in H1 and L1 containing this event, partitioned in 32 s segments. With a template bank of $2^{19}$ templates, a single spectrogram of the entire T duration of data requires cross-correlations for $64\times 2^{19}$. This is equivalent to 373 tera-floating point operations (teraFLOPs) in single precision (SP, f32), assuming that FFT is used and the data are sampled at 4096 Hz. This will be increased to 60 petaFLOPs when extending over $Q=161$ time slides (below). Calibration over 56 injections requires some 3.3 exaFLOPs to produce response curves in Figure 6.

The gain in evaluating correlations in the Fourier domain when array sizes are sufficiently large is well-known. A further gain is obtained by heterogeneous computing by offloading inverse FFTs to GPUs, realizing the speed-up of at least an order of magnitude (Figure 7). For matched filtering, this requires care in circumventing bandwidth limitations of the peripheral computer interface express (PCIe) between GPU(s) and the host CPU. The key involves the use of pre- and post-callback functions, with the latter limiting the GPU output to the tail of the matched filtering output, which may contain a candidate signal, discarding the bulk of the output [102].

Invariably, this approach involves some mixed-language programming, e.g., F90/C++ for the host and C99 for the GPU kernels. Open compute language (OpenCL) is a commonly used interfacing between the host and GPU [159,160]. Quite generally, the gain in the compute performance when offloading tasks to GPUs depends on the arithmetic intensity (e.g., [161]), expressing a ratio of demand on the compute versus memory resources. In practice, the latter includes limitations in the PCIe bandwidth between the GPU and host. The computed performance of the FFT exhibits a typical break as a function of the array size due to a loss in arithmetic intensity when the array size exceeds the limit of the local memory of the GPU’s compute units. Across the board, the performance drops well below the theoretical peak-compute performance, limited by memory bandwidth to global memory; access is relatively slow when a large number of clock cycles are involved.

Figure 7 shows variations in the arithmetic intensity of the OpenCL implementation clFFT on GPUs obtained in a benchmark. When array sizes are large, the results are about 7% of the theoretical peak-compute performance, showing the need for high-bandwidth memory (HBM) and, to scale acceleration, multiple GPU nodes. Circumventing the aforementioned PCIe bandwidth limitations is essential in the application of matched filtering.

6.2. Dynamical Load Balancing by Synaptic Parallel Processing

When analyzing a large number of data frames, high throughput can be realized by distributed computing on a platform of multi-GPU nodes organized over one or more local area networks (LANs). In the face of a heterogeneous platform consisting of nodes with varying compute performances, further optimization is obtained with dynamical load balancing by use of synaptic parallel processing [162]. This approach takes advantage of the embarrassingly parallel processing of a collection of independent tasks. Independent tasks can be specified in a list of instructions given by line items for dissemination on a platform of heterogeneous compute nodes organized in one or more local area networks (LANs).

Dynamical load-balanced processing is realized through synaptic computing, where individual nodes independently request line items from a shared task list with atomic access [162]. Access may be passive or through a traditional handshake protocol. This approach does not involve a central server, apart from a network-attached storage (NAS) device to collect output from all tasks onto a central location. In fact, this approach is the reverse of traditional approaches, where a server assigns tasks to clients.

7. Parameter Estimation in an Extended Foreground of (H1,L1)-Spectrograms

An event-timing analysis requires high-resolution parameter estimation by extracting probability density functions (PDFs) of observables, such as the start time $t_s$ and related parameters in (24). The expectation value and uncertainty can be inferred by calculating the mean and standard deviations from the data, which are needed for inferring the statistical significance of a candidate feature based on (11) and, independently, the consistency of observations in the independent detector channels (H1 and L1). Below, we discuss the PDFs obtained from H1 and L1, specifically the mean values from the merged (H1,L1) spectrograms. The independent H1 and L1 analysis is deferred to the following section.

For the descending chirp in gravitational radiation, a PDF($t_s,{\tau }_s$) is extracted, based on the parameter estimation (24) applied to the merged (H1,L1) spectrograms, using an indicator function $\chi (t_s,{\tau }_s,f_0,f_1)$. Our indicator function measures the correlation strength between H1 and L1 [81,109] as a function of all four parameters in (24). The evaluation involves heterogeneous computing and, to save memory, the output is limited to $\widehat{\chi }(t_s,{\tau }_s)$ at maxima over ($f_0,f_1)$.

7.1. PDFs from the Extended Foreground Analysis over Small Time Slides

A traditional foreground analysis with zero time slides between H1 and L1 data defines $\widehat{\chi }(t_s,{\tau }_s)$. A PDF($t_s,{\tau }_s)$ is obtained by repeating this estimate over a large number ($N_c=99$) of small time slides $\Delta t$ between H1 and L1 over the extended foreground

$$\begin{array}{c}\hfill \left|\Delta t\right|<\tau .\end{array}$$

(26)

Time slides $\Delta t$ are selected from a distribution of small and large time slides distributed symmetrically around zero, including a neighborhood of 0 (by 41 tiny time slides within the light travel time of 10 ms between H1 and L1).

The PDF($t_s)$ and PFF(${\tau }_s$) are derived from PDF($t_s,{\tau }_s)$ by tracing out ${\tau }_s$, respectively, $t_s$, i.e., the projection of a cluster defined by $\widehat{\chi }(t_s,{\tau }_s)$ in the $(t_s,{\tau }_s)$ plane of the parameter space over (26).

A candidate signal is identified by the cluster centered about $(t_s^{*},{\tau }_s^{*}$), defined by the maximum ${\widehat{\chi }}^{*}=\widehat{\chi }(t_s^{*},{\tau }_s^{*})$ of $\widehat{\chi }(t_s,{\tau }_s)$ over (26), as shown in Figure 8. (When estimating $(t_s^{*},{\tau }_s^{*}$), $(f_s,f_0)$ are the nuisance parameters projected out by taking maxima, as $\chi -$image analysis runs over all four parameters for each time slide $\Delta t$). The degree of clustering in the $(t_s,{\tau }_s)$ space is quantified by the number of elements in a box

$$\begin{array}{c}\hfill \left|t_s-t_s^{*}\right|<\tau ,\left|{\tau }_s-{\tau }_s^{*}\right|<\tau ,\end{array}$$

(27)

normalized to the total number $N_c$ of time slides (26). A 100% result defines an extremal cluster in the parameter space ($t_s,{\tau }_s)$.

For background, ensembles of peaks ${\widehat{\chi }}^{*}$ (and their degrees of clustering) are identified in each of the 1024 cells with a duration of $w=2$ s, covering an analysis over a total duration of $T=2024$ s using the aforementioned snippets of H1L1 data. Each cell has a total of $Q=161$ time slides over $[-32, 32]$ s, including (26). According to (26) and (27), up to $N_c$ can be clustered.

7.2. Clustering in Parameter Space

Absent of a gravitational-wave signal, $t_s$ and ${\tau }_s$ are uncorrelated and fail to cluster with a uniform statistic in $t_s$. In response to a gravitational-wave signal passing by, however, $t_s$ and ${\tau }_s$ cluster about a central value, which is representative of a correlation between H1 and L1 (Figure 8). The degree of clustering arising from such correlation strength is representative of the signal strength that, if present, will be persistent over the extended foreground (26). Crucially, extremal clustering provides a robust guard against false positives: weak signals or glitches occurring independently in H1 or L1 inevitably generate local maxima in $\widehat{\chi }$, but fail to produce extremal clustering.

Figure 9 evidences the uniqueness of extremal clustering in the extended emission feature, identifying

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

(28)

and satisfying the aforementioned causality condition (11). In particular, (11) is satisfied by PDF($t_s)$ (upon projecting out ${\tau }_s$) with

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

(29)

The start time $t_sϵG$ of the descending chirp in Figure 2 defines a trigger of GRB170817A, establishing the descending chirp as the emission from the central engine of GRB170817A in (3) [81]. In fact, Ref. (29) satisfies the more stringent causality condition in Butterfly matched filtering,

$$\begin{array}{c}\hfill C_1^{\prime }:t_m+\tau <t_s<t_m+t_g,\end{array}$$

(30)

where $\tau \simeq 0.5$ s refers to the duration of the templates in our Butterfly matched filtering.

The first PFA derives from $t_s$, satisfying $C_1$ by mere chance under the null hypothesis of a background of noise only with no gravitational wave passing by. According to (14), we have [118]

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

(31)

It follows that $C_1$ carries a PFA $p_1$ with an equivalent two-sided Gaussian-equivalent significance of $3.33\sigma$ and a false alarm rate (FAR) of about 1 per month. ($\mathrm{FAR}=\mathrm{PFA}/T$ (e.g., [163])). Box 4 summarizes the high-resolution observations on event timing in GRB170817-GRB170817A.

Box 4. High-resolution event timing in the extended foreground by the exascale HPC.

• Exascale HPC is achieved by synaptic parallel processing on a heterogeneous compute platform [162] in the mixed F90/C++/C99 using OpenCL and bash.
• Gavitational-wave observations in (3) provide event timing in high-resolution PDF($t_s$) by Butterfly matched filtering over an extended foreground over small time slides (10).
• The estimated $t_s$ satisfies the causality condition (11) and the more restricted condition (30).

8. Consistent Event Timing in Independent H1 and L1 Analyses

A further criterion for observation is the consistency in event timing in independent H1 and L1 analyses:

$$C_2:\left[t_s\right]=t_{s,H1}-t_{s,L1}$$

(32)

is zero within the resolution of our observation. This is evaluated by their respective PDFs$\left(t_s,{\tau }_s\right)$. The event timing analysis between H1 and L1 yields a PFA that is statistically independent of the previous PFA (14). Recall that the mean and differences of two estimates represent a rotation of the two-channel data over $\pi /4$ which corresponds to a unitary transformation of H1 and L1 data. The total PFA is obtained by factoring over these two PFAs.

We note that the joint PDF over $\left(t_{s,H1},t_{s,L1},{\tau }_{s,H1},{\tau }_{s,H1}\right)$ has an equivalent expression over $\left(t_s,{\tau }_s,\left[t_s\right],\left[{\tau }_s\right]\right)$, where $\left(t_s,{\tau }_s\right)$ is the mean and square brackets refer to the differences: $\left[t_s\right]=t_{s,H1}-t_{s,L1}$, $\left[{\tau }_s\right]={\tau }_{s,H1}-{\tau }_{s,L1}$. In the previous section, the mean is effective in the merged (H1,L1)-spectrograms.

With no gravitational wave passing by, H1 and L1 detectors are uncorrelated. In this case, the mean and difference in $\left(t_s,{\tau }_s\right)$ are statistically independent. The aforementioned PDF then factorizes as $p_A\left(t_s,{\tau }_s\right)p_B\left(\left[t_s\right],\left[{\tau }_s\right]\right)$. We are at liberty to project out differences, leaving $p_A\left(t_s,{\tau }_s\right)$ and (31), further projecting out ${\tau }_s$ based on (30). Projecting out $\left(t_s,{\tau }_s\right)$ reduces the same to $p_B\left(\left[t_s\right],\left[{\tau }_s\right]\right)$ independent of the context; that is, blind to GW170817-GRB170817 under the null hypothesis in (3).

Consistency is measured by cross-correlating the PDF($t_s$) obtained for H1 and L1. This gives rise to a PFA in the consistency of the time difference with zero, taking into account finite resolution in the two PDFs of $t_s$.

PDF($t_s$)s are obtained from H1 and L1 at high resolution, about an order of magnitude below $\tau$, by sampling over a large number of trials. (In the case at hand, the resolution remains above the light travel time between the two LIGO detectors). For this purpose, PDF($t_s$)s are generated by $nm$ trials: $n=10$ template seeds in Butterfly matched filtering with output branched into $m=32$ H1 and L1 spectrograms over the offsets $0,1,\cdots$, up to stride m.

Butterfly matched filtering uses a dense rather than orthogonal bank of templates, with each possibly providing an approximate match over the intermediate time scale of the phase coherence $\tau$ [101,102]. The output of data samples are typically covered by a multiplicity of neighboring templates. On average, time-stepping

$$\begin{array}{c}\hfill \Delta t_{ik}=t_{i_{k+1}}-t_{i_k}\end{array}$$

(33)

in $t_i$ is hereby below the sample time of $1/4096$ s of the H1L1 data, giving $\Delta t_{ik}$, which is typically 0 or $1/4096$ s, with STD ≃ $40\mathsf{\mu }$s. Therefore, on average, there are 7 samples in a stride of $m=32$, with a small STD of the order of $\mathcal{O}\left({10}^{-2}\right)$ Hz [109]. The high resolution is derived from offsets $m_0=1,2,\cdots ,m$, effectively, by subsampling.

Figure 10 and Figure 11 show the results of independent H1 and L1 analyses by PDF($t_s$) for each one. For the total duration $T=2048$ s of H1L1 data, the local maximum of $\widehat{\chi }\left(t_s,{\tau }_s\right)$ shows a correlation of ${\rho }_{\left[t_s\right],\left[{\tau }_s\right]}=1.29\times {10}^{-4}$, where $\left[t_s\right]=t_{s,H1}^{*}-t_{s,L1}^{*}$ and $\left[{\tau }_s\right]={\tau }_{s,H1}^{*}-{\tau }_{s,L1}^{*}$. This is expected in independent H1 and L1 analyses [109]. PDFs($t_s$) of H1 and L1 are maximally correlated at a time slide that is consistent with zero within the resolution $\delta t_{res}=0.025$s of our binning. This maximum is global; it is maximal among all the maxima identified in each of the 64 segments of 32 s, covering the entire analysis over $T=2048$s. By (32),

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

(34)

is the $p$-value of consistency with two-sided Gaussian-equivalent significance of $4.06\sigma$ and the FAR of ∼1 per 1.4 years. The factor of 4 in (34) is a combination of 2 coefficients. The first is a factor of 2, which originates from taking the absolute value in $\left|\left[t_s\right]\right|<\delta t_{res}$ at the peak of the cross-correlation between PDF$\left(t_{s,H1}\right)$ and PDF$\left(t_{s,L1}\right)$. Additionally, the uniform prior in astrophysical timing $t_{H1}$ and $t_{L1}$ leads to a top-hat distribution in $\left[t_s\right]$, giving another factor of 2.

In the mean of the PDF($t_s$)s of H1 and L1, the delay time in the gravitational collapse to a black hole satisfies

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

(35)

consistent with (29), in the mean from the merged (H1,L1)-spectrograms. The delay times (35) and (29) identify the trigger of GRB170817A by the associated descending chirp in gravitational radiation. It is notably consistent with the lifetime $t_w=0.{98}_{-0.26}^{+0.31}$ s of the progenitor HMNS derived from (neutron-rich) outflow propagation times and the mass of blue ejecta [129]; its survival time is based on neutrino cooling and internal dissipation of rotational energy (e.g., [119]). Maximal clustering (35) comes with PDF$\left({\tau }_s\right)$. Mean values of ${\tau }_s$ in H1 and L1 are, respectively, ${\tau }_s=\left(3.1\pm 0.1\right)$s and ${\tau }_s=\left(2.9\pm 0.2\right)$s. ${\tau }_s=\left(3.1\pm 0.1\right)$s (H1) and ${\tau }_s=\left(2.9\pm 0.2\right)$s (L1) associated with the cluster (35) satisfy the mean

$$\begin{array}{c}\hfill {\tau }_s=\left(3.00\pm 0.09\right)\mathrm{s}.\end{array}$$

(36)

The result (36) is consistent with ${\tau }_s$, derived from (29) in the merged (H1,L1)-spectrograms; (29) and (36) are representative of the lifetime of the black hole spin and provide constraints on $T_{90}^{8-70\mathrm{keV}}$ of GRB170817A [90,109]. Box 5 summarizes the two independent analyses.

Box 5. Statistically independent PFAs of the trigger time $t_s$ of GRB170817A.

• Independent event timings derive from the mean and difference in $t_s$, equivalent to a unitary transformation of H1 and L1 data, here obtained from merged and individual H1 and L1 spectrograms.
• PFA${}_1$ (14) derives from causality (11) satisfied in merged spectrograms. A second PFA${}_2$ (34) derives from consistency in individual H1 and L1 spectrograms.
• A total PFA with an equivalent Z-score of 5.5 is factored over the ordinary probability (14) and the p-value (34).

9. Conclusions and Outlook

The EM $\nu$ observations of SN1987A demonstrate the power of multi-messenger observations in identifying the formation of an NS that would otherwise be hidden by the opaque stellar envelope. In the absence of a rigorous identification of the nature of the remnant [164], however, the true trigger of the explosion eludes this neutrino observation. It is unknown whether the NS underwent further collapse to become a BH.

9.1. Principal Results

The EM-GW observations of GRB170817A provide evidence for the progenitor system of a short GRB: the DNS system GW170817. In this case, there is an additional observation of a trigger in the gap of 1.7 s between GW170817 and GRB170817A. It signals the start time $t_s$ of a descending chirp in gravitational radiation characteristics of the spin-down of a compact object powering the central engine of GRB170817A, with extended emissions covering the entire duration of this burst (Box 1). This is the first detection of a non-merger gravitational wave signal, discovered in high-resolution spectrograms (Box 2). Note that the finite delay time (35) in gravitational collapse to a Kerr black hole (Box 3), in agreement with independent estimates of the lifetime of the HMNS (Figure 2), rules out a BH-NS merger progenitor to GRB170817A [165].

PFA${}_1$ (31) and PFA${}_2$ (34) pertain to the mean and difference in $t_s$ based on the merged (H1,L1) spectrograms and, respectively, individual H1 and L1 spectrograms obtained from event timing constraints $C_{1,2}$, evaluated by statistically independent PDFs($t_s)$ therein. The PDFs($t_s)$ of $t_s$ are extracted by means of HPC, accelerated by software, hardware, and distributed computing (Box 4). A Z-score of 5.49 is derived from the joint PFA $p_1\times p_2=4.05\times {10}^{-8}$, and is on par with the GW170817-GRB170817A association [6] (Box 5).

Table 1. Event timing of ${\mathcal{E}}_{GW}$ (5) in (3) with type I error statistics, independently in causality ($C_1$: $p_1$) and consistency ($C_2$: $p_2$) in independent H1 and L1 analyses over a duration T. Central values and uncertainties refer to the mean and standard deviations of PDFs. (Reprinted from [109]).

H1	L1	H1	L1	H1,L1	H1,L1	Merged (H1,L1)	Merged (H1,L1)
$t_s-t_m$ [s]	$t_s-t_m$ [s]	${\tau }_s$ [s]	${\tau }_s$ [s]	$t_s-t_m$ [s]	${\tau }_s$ [s]	$t_s-t_m$ [s]	${\tau }_s$ [s]
$0.9130\pm 0.1366$	$0.9234\pm 0.1122$	$3.1\pm 0.1$	$2.9\pm 0.2$	$0.92\pm 0.08$	$3.00\pm 0.09$	$0.86\pm 0.10$	$2.91\pm 0.17$
T [s]	$p_1$	FAR${}_1^{-1}$	$p_2$	FAR${}_2^{-1}$	$p_1\times p_2$	FAR${}^{-1}$
2048	$1.7/2048$	1 month	$0.1/2048$	1.4 yr	$4.05\times {10}^{-8}$	1.6 kyr
204,800	1.7/204,800	782 yr	$0.1/2048$	1.4 yr	$4.05\times {10}^{-10}$	>$160$ kyr

summarizes these PFAs.

The first direct observation of the spin-down of a Kerr BH is of principal significance to extreme transient events, by the universality of black holes in modern astronomy and the unique power of Kerr black holes compared to neutron stars. The high confidence level of this observation (Table 1) and the agreement in event timing in EM-GW (Box 1 and Box 3) suggest that this process may be observed in upcoming LVK observational runs.

Figure 12 presents a chronicle of (3) summarizing the observational results.

9.2. Outlook on Upcoming LVK Observations

The present results have potentially significant implications for upcoming observations of NS-NS mergers. These mergers inevitably have masses that are different from GW170817, resulting in commensurably different lifetimes of the initial HMNS produced in the coalescence. Observing multiple DNS mergers (sufficiently nearby) offers the opportunity to measure the correlation between the total system mass and the lifetime of the initial HMNS (reviewed for O4 in [149]), providing a potentially powerful constraint of the equation of state of neutron star matter.

In [81], the descending chirp is observed to agree remarkably well with quadrupole gravitational radiation emitted at a radius of $K\simeq 4$ times the ISCO around the Kerr BH formed in the gravitational collapse of the initial HMNS at the time (35). Since black holes do not retain the memory of their progenitors, except for the total mass and angular momentum, these findings are particularly interesting in the context of type Ib/c supernovae that give rise to black holes. This is also relevant in relation to their connection (albeit with a small branching ratio) to normal long GRBs [67] and superluminous supernovae [166,167].

LVK O4-5 observations are more favorable for detecting joint EM-GW signals from SN Ib/c, considering their relatively larger mass compared to the remnants of DNS mergers. A high black hole mass provides an advantage in terms of the potential for enhanced energy emissions at relatively low frequencies. These low frequencies align with or may overlap with the bandwidth of maximal sensitivity of LVK detectors. While rare, a failed GRB-SNe may, nevertheless, be luminous in gravitational radiation, and more frequent than DNS mergers, which may be tested with upcoming LVK probes in SN Ib/c events within the Local Universe.