Search for Gravitational Waves from a Long-lived Remnant of the Binary Neutron Star Merger GW170817

STAMP (Thrane et al. 2011) is an unmodeled search pipeline that is designed to detect gravitational wave transients. Its basic unit is a spectrogram made from cross-correlated data between two detectors. Narrowband transient gravitational waves produce tracks of excess power within these spectrograms and can be detected by pattern-recognition algorithms. Each spectrogram pixel is normalized with the noise to obtain a signal-to-noise ratio (S/N) for each pixel.

STAMP was used in the first GW170817 post-merger search (Abbott et al. 2017g) in a configuration with 500 s long spectrograms. To increase sensitivity to longer GW signals, here we use spectrogram maps of 15,000 s length. The search is split into two frequency bands from 30 to 2000 and 2000 to 4000 Hz. The former uses pixels of $100\,{\rm{s}}\times 1\,\mathrm{Hz}$, while the latter uses shorter-duration pixels of $50\,{\rm{s}}\times 1\,\mathrm{Hz}$ to limit S/N loss due to the Earth's rotation changing the phase difference between detectors.

We then use Stochtrack (Thrane & Coughlin 2013), which is a seedless clustering algorithm, to identify significant clusters of pixels within these maps. The algorithm uses one million quadratic Bézier curves as templates for each map and the loudest cluster is picked for each map. More details about the pixel size choice, the detection statistic and the search results are given in Appendix A.1.

The on-source data window is from just after the time of the merger to the end of O2 (1187008942–1187733618). To measure the background and estimate the significance of the clusters that we have found, we run the algorithm on time-shifted data from June 24th to just before the merger.

Hidden Markov model (HMM) tracking provides a computationally efficient strategy to detect and estimate a quasimonochromatic GW signal with unknown frequency evolution and stochastic timing noise (Suvorova et al. 2016; Sun et al. 2018). It was applied to data from the first aLIGO observing run to search for CWs from the low-mass X-ray binary Scorpius X-1 (Abbott et al. 2017f). The revision of the algorithm in Sun et al. (2018) is also well-suited to searching for a long-transient signal from a BNS merger remnant, if the spin-down timescale is in the range ${10}^{2}\,{\rm{s}}\lesssim \tau \lesssim {10}^{4}\,{\rm{s}}$.

A HMM is an automaton that is based on a Markov chain (a stochastic process transitioning between discrete states at discrete times), which is composed of a hidden (unmeasurable) state variable and a measurement variable. A HMM is memoryless; i.e., the hidden state at time ${t}_{n+1}$ only depends on the state at time t_n, with a certain transition probability. The most probable sequence of hidden states given the observations is computed by the classic Viterbi algorithm (Viterbi 1967). Further details on the probabilistic model can be found in Appendix A.2.

In this analysis, we track the GW signal frequency as the hidden variable, whose discrete states are mapped one-to-one to the frequency bins in the output of a frequency-domain estimator computed over an interval of length ${T}_{\mathrm{drift}}$. We aim to search for signals with ${10}^{2}\,{\rm{s}}\lesssim \tau \lesssim {10}^{4}\,{\rm{s}}$, such that the first time derivative ${\dot{f}}_{\mathrm{gw}}$ of the signal frequency ${f}_{\mathrm{gw}}$ satisfies ${\dot{f}}_{\mathrm{gw}}\,\approx {f}_{\mathrm{gw}}/\tau \lesssim 1\,\mathrm{Hz}\,{{\rm{s}}}^{-1}$, given ${T}_{\mathrm{drift}}=1\,{\rm{s}}$ and a frequency bin width of ${\rm{\Delta }}f=1\,\mathrm{Hz}$. The motion of the Earth with respect to the solar system barycenter (SSB) can be neglected during a ${T}_{\mathrm{drift}}$ interval. Hence, we use a running-mean normalized power in SFTs with length ${T}_{\mathrm{SFT}}={T}_{\mathrm{drift}}=1\,{\rm{s}}$ as the estimator to calculate the HMM emission probability.

We analyze 9688 s of data (GPS times 1187008882–1187018570) in a 100–2000 Hz frequency band with multiple configurations optimized for different τ. We do not analyze longer data stretches because: (i) several intervals in the data after GPS time 1187018570 are not in analyzable science mode, and (ii) signals with ${10}^{2}\,{\rm{s}}\lesssim \tau \lesssim {10}^{4}\,{\rm{s}}$ drop below the algorithm's sensitivity limit after $\sim {10}^{4}\,{\rm{s}}$. Observing for longer merely accumulates noise without improving S/N. The 9688 SFTs are Hann-windowed. The detection statistic ${ \mathcal P }$ is defined in Equation (17). The methodology and analysis is fully described in Sun & Melatos (2018).

The Adaptive Transient Hough search method is described in detail in Oliver et al. (2019). It follows a semi-coherent strategy similar to the SkyHough (Krishnan et al. 2004; Sintes & Krishnan 2007; Aasi et al. 2014) all-sky CW searches, but adapted to rapid-spindown transient signals.

We start from data in the form of Hann-windowed SFTs with lengths of [1, 2, 4, 6, 8] s, covering one day after merger (GPS times 1187008882–1187095282). These are digitized by setting a threshold of 1.6 on their normalized power, as first derived by Krishnan et al. (2004), replacing each SFT by a collection of zeros and ones called a peak-gram. For each point in parameter space, the Hough number count is the weighted sum of the peak-grams across a template track accounting for Doppler shift and the spindown of the source. The use of weights minimizes the influence of time-varying detector antenna patterns and noise levels (Sintes & Krishnan 2007). For this post-merger search, it also accounts for the amplitude modulation related to the transient nature of the signal.

The search parameter space for the model from Section 2 covers a band of 500–2000 Hz in starting frequencies ${f}_{\mathrm{gw}0}$, braking indices of $2.5\leqslant n\leqslant 7$ and spindown timescales of ${10}^{2}\leqslant \tau \leqslant {10}^{5}\,{\rm{s}}$. The search runs over 16,042 subgroups, each containing a range of 150 Hz in f₀, 0.25 in n and 1000 s in τ. Each subgroup is analyzed with the longest possible SFTs according to the criterion (Oliver et al. 2019)

$$\begin{eqnarray}&&{T}_{\mathrm{SFT}}\leqslant \displaystyle \frac{\sqrt{(n-1)\tau }}{\sqrt{{f}_{\mathrm{gw}}}},\end{eqnarray} \tag{ 5 }$$

and for each template the observation time is selected as ${T}_{\mathrm{obs}}=\min (4\,\tau ,24\,\mathrm{hr})$. Over the whole template bank, the search uses data from 187 to 2000 Hz.

Each template is ranked based on the deviation of its weighted number count from the theoretical expectation for Gaussian noise (the critical ratio) as described in Appendix A.3. The detection threshold corresponds to a two-detector 5σ false-alarm probability for the entire template bank. A per-detector critical ratio threshold was also set to check the consistency of a signal between H1 and L1.

The FrequencyHough is a pattern-recognition technique that was originally developed to search for CWs by mapping points in time-frequency space of the detector to lines in frequency-spindown space (Antonucci et al. 2008; Astone et al. 2014). This only works if the signal frequency varies in time very slowly. Miller et al. (2018) have generalized the FrequencyHough for post-merger signals, where we expect much higher spindowns.

The search starts at a time offset ${\rm{\Delta }}t={t}_{\mathrm{start}}-{t}_{{\rm{c}}}$ after coalescence time ${t}_{{\rm{c}}}$, so that the waveform model is interpreted with starting frequency ${f}_{\mathrm{start}}={f}_{\mathrm{gw}}(t={\rm{\Delta }}t)$ taking the place of ${f}_{\mathrm{gw}0}$ in Equation (3). In this way, assuming that the NS has already spun down before ${t}_{\mathrm{start}}$ following some arbitrary track, we would be able to probe higher initial frequencies and spindowns through a less challenging parameter space during the search window. Furthermore, the source parameters $(n,{f}_{\mathrm{start}},\tau$) are transformed to new coordinates, such that in the new space the behavior of the signal is linear. See Appendix A.4 for the transformation relations.

We search across the parameter space with a fine, nonuniform grid: for each braking index n, we do a Hough transform and then record the most significant candidates over the parameter range of the resulting map. This is done separately on the data from each detector. We then check the candidates for coincidence between detectors according to their Euclidean distance in parameter space.

The search is run in three configurations using varying ${T}_{\mathrm{SFT}}=2,4,8\,{\rm{s}}$, covering different observing times, starting ${\rm{\Delta }}t=1\mbox{--}7$ hr after merger. It covers $n=[2.5,7]$, ${f}_{\mathrm{start}}\,=[500,2000]$ Hz and $\tau =[10,{10}^{5}]\,{\rm{s}}$, analyzing detector data from 50 to 2000 Hz.

Candidates are also ranked by critical ratio (deviation from the theoretical expectation for Gaussian noise) in this analysis. Most can be vetoed by the coincidence step or by considering detector noise properties. A follow-up procedure for surviving candidates is also described in Appendix A.4.

All of the four search methods either found no significant candidates in the aLIGO data after GW170817; or those that were found, were clearly vetoed as instrumental artifacts.

For the unmodeled STAMP search, the loudest triggers in the low- and high-frequency bands have S/Ns of 3.18 and 3.07, respectively. The time-shifted backgrounds only just start to drop off near these S/Ns, so that they correspond to false-alarm probabilities ${p}_{\mathrm{FA}}$ of 0.81 and 0.80, which are completely consistent with noise. For reference, ${p}_{\mathrm{FA}}=0.05$ would only have been reached for S/Ns of $\gtrsim 4.9$ and $\gtrsim 3.5$ for these low- and high-frequency background distributions, respectively. (See Figure 4 in Appendix A.1.)

For HMM, the loudest trigger has a detection statistic ${ \mathcal P }=2.6749$ (as defined in Equation (17)), which corresponds to a false-alarm probability of 0.01, right below the threshold set beforehand as significant enough for further study. The trigger is found with observing time ${T}_{\mathrm{obs}}=200\,{\rm{s}}$ starting from $t={t}_{{\rm{c}}}$. Monte-Carlo simulations show that for the signals that this setup is sensitive to, a higher ${ \mathcal P }$ should be obtained with longer ${T}_{\mathrm{obs}}$. Follow-up analysis of the trigger with $300\,{\rm{s}}\,\leqslant {T}_{\mathrm{obs}}\leqslant 1000\,{\rm{s}}$ confirms that it does not follow this expectation; hence, it is discarded as spurious.

ATrHough found 51 initial candidates over the covered part of $(n,{f}_{\mathrm{gw}0},\tau )$ parameter space. All of these were excluded with the follow-up procedure described in Appendix A.3 because they are inconsistent with the expected spindown model and are more likely to have been caused by monochromatic detector artifacts (lines) contaminating the search templates.

The FreqHough search returned 521 candidates over the covered part of $(n,{f}_{\mathrm{start}},\tau )$ parameter space. We vetoed 10 of them because they were within frequency bands contaminated with known noise lines (Covas et al. 2018). A total of 510 of the remaining candidates had much higher ($\gt 4$ times) critical ratios in H1 than in L1, which is inconsistent with true astrophysical signals when considering the relative sensitivities, duty factors and antenna patterns. There was one remaining candidate, with a critical ratio of 5.21 in H1 and 4.88 in L1, which was followed up and excluded with the procedure described in Appendix A.4.

Starting from this non-detection result, we use simulated signals according to Equation (3) to quantify the sensitivity of each analysis given the dataset around the time of GW170817 and its known sky location. The sets of injected parameters are different for each pipeline, while there are also some differences in procedure: STAMP performs injections on the same data as the main search but with a non-physical time shift between the detectors (as in Abbott et al. 2017g, 2018a); HMM injects signals into the original set of SFTs but with randomly permuted timestamps; and the other two pipelines inject signals into exactly the same data as analyzed in the main search. HMM and ATrHough perform all injections starting at merger time ${t}_{{\rm{c}}}$, with ${f}_{\mathrm{gw}0}$ in Equation (3) interpreted as the frequency at ${t}_{{\rm{c}}}$, while injections for FreqHough are done at ${\rm{\Delta }}t=1,2$ or 5 hr after ${t}_{{\rm{c}}}$, which are chosen as representative starting times for each search configuration, and ${f}_{\mathrm{gw}0}$ is correspondingly set at ${t}_{{\rm{c}}}+{\rm{\Delta }}t$. Similarly, STAMP treats ${f}_{\mathrm{gw}0}$ as the starting frequency of each injection, which have ${\rm{\Delta }}t$ distributed through the whole search range, yielding a time-averaged sensitivity. In the following, we use ${f}_{\mathrm{start}}$ to refer to any of these choices.

These differences in injection procedure and the different choices of detection threshold mean that any comparison of the following results does not correspond to a representative evaluation of general pipeline performance but is solely in the interest of estimating how much sensitivity is missing for a GW170817-like post-merger detection that is based on the specific configurations that are used in the present search.

We focus here on results for a braking index of n = 5, as expected for spin-down dominated by GW emission from a static quadrupole deformation. The signal amplitude h₀ (as given in Equation (4)) is degenerate between the ellipticity , moment of inertia ${I}_{{zz}}$ and distance d. We choose a fiducial value of ${I}_{{zz}}\,=100\,{M}_{\odot }^{3}\,{G}^{2}/{c}^{4}\approx 4.34\times {10}^{38}\,\mathrm{kg}\,{{\rm{m}}}^{2}$, which is consistent with EoS yielding a supramassive or stable remnant: the high mass and assumed rapid rotation can increase the moment of inertia by more than a factor of 3 compared to a nonrotating NS of $1.4\,{M}_{\odot }$. In addition, EoS compatible with the high remnant mass favor larger moments of inertia already at lower mass. For a given set of model parameters $\{n=5,{f}_{\mathrm{start}}\,,\tau \}$, we consider the maximum allowed by the initial rotational energy budget (Sarin et al. 2018): the total emitted GW energy as $t\to \infty$,

$$\begin{eqnarray}&&{E}_{\mathrm{gw}}=-{\int }_{t={t}_{\mathrm{start}}}^{\infty }{dt}\,\displaystyle \frac{32G}{5{c}^{5}}\,{I}_{{zz}}^{2}\,{\epsilon }^{2}\,{{\rm{\Omega }}}^{6}(t),\end{eqnarray} \tag{ 6 }$$

must not exceed the remnant's initial rotational energy ${E}_{\mathrm{rot}}\,=0.5\,{I}_{{zz}}\,{f}_{\mathrm{start}}^{2}\,{\pi }^{2}$.

Given each pipeline's detection threshold, we can rescale the amplitude of simulated signals until 90% of them are recovered above threshold, while randomizing over nuisance parameters (polarization angle and initial phase of the signal; also source inclination ι and signal start time for some of the pipelines). We can then either interpret this amplitude scaling as a need to lower the distance of simulated sources; i.e., estimating the sensitive distance ${d}^{90 \% }$ of the search. Alternatively, we can fix the true distance to the source of GW170817 to obtain an energy upper limit ${E}_{\mathrm{gw}}^{90 \% }$: for a remnant NS with the given parameters, we interpret the square of the amplitude scaling as the factor by which the energy output needs to be greater than the expected amount to produce signals that we can recover at 90% confidence. Both interpretations are shown in Figure 2, with coverage of the injection sets illustrated in Figure 3. The full results are listed in the appendix in Tables 2–5.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The highest sensitivities are achieved at low ${f}_{\mathrm{start}}$ and for rapid spindown (low τ). This is mostly due to the energy budget constraint enforced on : in principle, higher ${f}_{\mathrm{start}}$ yield higher initial amplitudes and longer τ allow accumulation of S/N over longer observation times (${\rm{S}}/{\rm{N}}\propto {\tau }^{1/2}$ for n = 5 and once ${T}_{\mathrm{obs}}\gt \tau$). However, due to the energy budget constraint, must be lower in this region of parameter space and hence actual detectability is reduced (${\rm{S}}/{\rm{N}}\propto \epsilon \propto {\tau }^{-1}{f}_{\mathrm{start}}^{-6}$).

The four pipelines perform differently across τ regimes: the unmodeled STAMP and HMM are most sensitive at the shortest $\tau ={10}^{2}$ s, but lose up to an order of magnitude in ${d}^{90 \% }$ when going to $\tau ={10}^{4}\,{\rm{s}}$. Meanwhile, the model-based semi-coherent ATrHough and FreqHough have focused on longer τ of $4\times {10}^{2}\,{\rm{s}}$ to $3\times {10}^{4}\,{\rm{s}}$, with only up to a factor of 2 loss in ${d}^{90 \% }$ for the longest τ at fixed ${f}_{\mathrm{start}}$. See Figures 6–9 in the appendix for sensitivity estimates over each pipeline's full injection set.

Given that this parameter dependence is shaped by the ${E}_{\mathrm{gw}}={E}_{\mathrm{rot}}$ constraint and is also influenced by some practical tradeoffs in pipeline configuration, in this paper we do not attempt to provide a general evaluation of pipeline performance on fully equivalent injection sets or for generic GW signals. Such a comparison would require a detailed mock data challenge similar to Messenger et al. (2015) and Walsh et al. (2016). Instead, Figure 2 shows results from each pipeline for the parts of parameter space where it achieved its highest sensitivity.

In summary, in no part of the n = 5 parameter space covered by the four search ranges and injection sets do we reach 90% sensitive distances of 1 Mpc or further. This corresponds to a lowest 90% detectable energy of ${E}_{\mathrm{gw}}\lesssim 8\,{M}_{\odot }{c}^{2}$ at ${f}_{\mathrm{start}}\,=500\,\mathrm{Hz}$ and $\tau =100\,{\rm{s}}$. At higher ${f}_{\mathrm{start}}$, the sensitive distances for any τ are even lower due to the energy constraint. Note again that this covers power-law spindown signals, both starting immediately at coalescence time ${t}_{{\rm{c}}}$ and signals starting with some time delay, with a delay time of 1–7 hr for FreqHough and any possible delay time until the end of O2 for STAMP.

At the shortest τ, the parameter space covered here overlaps with the magnetar injections in the shorter-duration search of Abbott et al. (2017g),¹⁷² although the results in that paper were quoted as recoverable at 50% confidence, and hence are more optimistic than the new results at 90%. For example, at ${f}_{\mathrm{start}}=1000\,\mathrm{Hz}$ and $\tau =100\,{\rm{s}}$, the STAMP analysis in the previous paper found ${E}_{\mathrm{gw}}^{50 \% }\approx 24\,{M}_{\odot }{c}^{2}$ while the new STAMP and HMM analyses presented here obtain ${E}_{\mathrm{gw}}^{90 \% }\approx 100\,{M}_{\odot }{c}^{2}$ at these parameters. For the pipelines in this paper, amplitudes for detectability at 50% confidence are typically lower by a factor of 2–4 than those at 90% confidence. Although these lower thresholds would push the best ${d}^{50 \% }$ limits up to a few Mpc, this would not change the conclusion that any GWs from a long-lived remnant of GW170817 at 40 Mpc would be undetectable.

Spectrogram pixel sizes. The low-frequency band from 30 to 2000 Hz uses pixels of $100\,{\rm{s}}\times 1\,\mathrm{Hz}$, while the high-frequency band from 2000 to 4000 Hz uses pixels of smaller durations of $50\,{\rm{s}}\times 1\,\mathrm{Hz}$. Smaller pixels at higher frequency are necessary to account for the rotation of the Earth, which causes the GW phase difference between detectors to change with time. There is a loss of S/N when the pixel durations are too large, which increases with frequency. Therefore, the durations are chosen to limit the maximum possible S/N loss in a pixel (at the highest frequencies) from this effect to about 10% (Thrane et al. 2015).

Detection statistic. Each spectrogram in the STAMP search (Thrane et al. 2011; Thrane & Coughlin 2013; Thrane et al. 2015) is analyzed with many randomly chosen quadratic Bézier curves. The S/N of each track ${\rho }_{{\rm{\Gamma }}}$ is a weighted sum of the S/N of the pixels covered by the track. The quantity ${\rho }_{{\rm{\Gamma }}}$ also serves as the detection statistic and is calculated as:

$$\begin{eqnarray}&&{\rho }_{{\rm{\Gamma }}}=\displaystyle \frac{1}{{N}^{3/4}}\displaystyle \sum _{i}{\rho }_{i},\end{eqnarray} \tag{ 7 }$$

where i runs over all the pixels in a track and N is the total number of pixels in it. These are then ranked and the track with largest ${\rho }_{{\rm{\Gamma }}}$ is picked as the trigger for a map. This is done for both the main on-source search and for the background estimation over time-shifted data.

Background triggers and loudest events. Figure 4 shows the distribution of false-alarm probabilities ${p}_{\mathrm{FA}}$ for the S/Ns of triggers collected in background data, for both high- and low-frequency spectrograms. The loudest on-source event in each frequency range is also shown.

Zoom In Zoom Out Reset image size

A general description of the HMM method is given by Suvorova et al. (2016) and Sun et al. (2018). The following summary is intended to clarify the configuration used for the search presented in this paper.

Probabilistic model. A Markov chain is a stochastic process transitioning between discrete states at discrete times $\{{t}_{0},\,\cdots ,\,{t}_{{N}_{T}}\}$. A HMM is an automaton based on a Markov chain, composed of the hidden (unmeasurable) state variable $q(t)\in \{{q}_{1},\,\cdots ,\,{q}_{{N}_{Q}}\}$ and the measurement variable $o(t)\in \{{o}_{1},\,\cdots ,\,{o}_{{N}_{O}}\}$. A HMM is memoryless; i.e., the hidden state at time ${t}_{n+1}$ only depends on the state at time t_n, with transition probability

$$\begin{eqnarray}&&{A}_{{q}_{j}{q}_{i}}=P[q({t}_{n+1})={q}_{j}| q({t}_{n})={q}_{i}].\end{eqnarray} \tag{ 8 }$$

The hidden state q_i is in observed state o_j at time t_n with emission probability

$$\begin{eqnarray}&&{L}_{{o}_{j}{q}_{i}}=P[o({t}_{n})={o}_{j}| q({t}_{n})={q}_{i}].\end{eqnarray} \tag{ 9 }$$

Given the prior defined by

$$\begin{eqnarray}&&{{\rm{\Pi }}}_{{q}_{i}}=P[q({t}_{0})={q}_{i}],\end{eqnarray} \tag{ 10 }$$

the probability that the hidden state path $Q=\{q({t}_{0}),\,\cdots ,\,q({t}_{{N}_{T}})\}$ gives rise to the observed sequence $O=\{o({t}_{0}),\,\cdots ,\,o({t}_{{N}_{T}})\}$ equals

$$\begin{eqnarray}&&P(Q| O)={L}_{o({t}_{{N}_{T}})q({t}_{{N}_{T}})}{A}_{q({t}_{{N}_{T}})q({t}_{{N}_{T}-1})}\cdots {L}_{o({t}_{1})q({t}_{1})}{A}_{q({t}_{1})q({t}_{0})}{{\rm{\Pi }}}_{q({t}_{0})}.\end{eqnarray} \tag{ 11 }$$

The most probable path

$$\begin{eqnarray}&&{Q}^{* }(O)=\mathrm{argmax}P(Q| O),\end{eqnarray} \tag{ 12 }$$

maximizes $P(Q| O)$ and gives the best estimate of q(t) over the total observation, where $\mathrm{argmax}(\cdots )$ returns the argument that maximizes the function $(\cdots )$. We use the classic Viterbi algorithm (Viterbi 1967) to efficiently solve the HMM and compute ${Q}^{* }(O)$.

Search setup. In this analysis, we track $q(t)={f}_{\mathrm{gw}}(t)$, where ${f}_{\mathrm{gw}}(t)$ is the GW signal frequency at time t. The discrete hidden states are mapped one-to-one to the frequency bins in the output of a frequency-domain estimator G(f) computed over an interval of length ${T}_{\mathrm{drift}}$, with bin size ${\rm{\Delta }}f$. We aim to search for signals with ${10}^{2}\,{\rm{s}}\lesssim \tau \lesssim {10}^{4}\,{\rm{s}}$, corresponding to ${\dot{f}}_{\mathrm{gw}}\lesssim 1\,\mathrm{Hz}\,{{\rm{s}}}^{-1}$. Therefore, we choose ${T}_{\mathrm{drift}}=1\,{\rm{s}}$ (i.e., ${\rm{\Delta }}f=1\,\mathrm{Hz}$) to satisfy

$$\begin{eqnarray}&&\left|{\int }_{t}^{t+{T}_{\mathrm{drift}}}{dt}^{\prime} {\dot{f}}_{\mathrm{gw}}(t^{\prime} )\right|\leqslant {\rm{\Delta }}f\end{eqnarray} \tag{ 13 }$$

for $0\leqslant t\leqslant {T}_{\mathrm{obs}}$, where ${T}_{\mathrm{obs}}$ is the total observing time. The motion of the Earth with respect to the SSB can be neglected during the interval [$t,t+{T}_{\mathrm{drift}}$]. Hence, the emission probability ${L}_{o(t){q}_{i}}=P[o(t)| {f}_{i}\leqslant {f}_{\mathrm{gw}}(t)\leqslant {f}_{i}+{\rm{\Delta }}f]\propto \exp [G({f}_{i})]$ is calculated from the running-mean (window width 3 Hz) normalized power in SFTs with length ${T}_{\mathrm{SFT}}={T}_{\mathrm{drift}}=1\,{\rm{s}}$ as the estimator G(f). We write

$$\begin{eqnarray}&&G({f}_{i})=\displaystyle \sum _{X}{\tilde{y}}_{i}^{X}{\tilde{y}}_{i}^{X* },\end{eqnarray} \tag{ 14 }$$

where i indexes the frequency bins of the normalized SFT $\tilde{y}$, X indexes the detector, and the repeated index i on the right-hand side does not imply summation. We assume that the auto-correlation timescale of timing noise is much longer than ${T}_{\mathrm{drift}}$, and hence adopt the transition probabilities

$$\begin{eqnarray}&&{A}_{{q}_{i-1}{q}_{i}}={A}_{{q}_{i}{q}_{i}}=\displaystyle \frac{1}{2},\end{eqnarray} \tag{ 15 }$$

with all other entries being zero. Since we have no independent knowledge of ${f}_{\mathrm{gw}}$, we choose a uniform prior, viz

$$\begin{eqnarray}&&{{\rm{\Pi }}}_{{q}_{i}}={N}_{Q}^{-1}.\end{eqnarray} \tag{ 16 }$$

We define a detection statistic ${ \mathcal P }$, given by

$$\begin{eqnarray}&&{ \mathcal P }=\displaystyle \frac{1}{{N}_{T}+1}\displaystyle \sum _{n=0}^{{N}_{T}}G[{f}_{i({t}_{n})}],\end{eqnarray} \tag{ 17 }$$

where the integer $i({t}_{n})$ indexes the SFT frequency bin corresponding to ${q}^{* }({t}_{n})$ on the optimal path Q* (${t}_{0}\leqslant {t}_{n}\leqslant {t}_{{N}_{T}}$).

The strain amplitude h₀ in (4) decreases significantly for $t\gg \tau$. Hence, the instant S/N decreases for ${T}_{\mathrm{obs}}\gtrsim \tau$. Monte-Carlo simulations show that choosing ${T}_{\mathrm{obs}}\sim \tau$ yields the best sensitivities for signals with h₀ (Equation (4)) near the detection limit for the waveform in Equation (3).

The initial spin-down rate $| {\dot{f}}_{\mathrm{gw}0}|$ of a signal with $\tau \lesssim {10}^{3}\,{\rm{s}}$ can be too high (i.e., $| {\dot{f}}_{\mathrm{gw}0}| \gt 1\,\mathrm{Hz}\,{{\rm{s}}}^{-1}$) for Equation (13) to be satisfied with ${T}_{\mathrm{drift}}=1\,{\rm{s}}$. We start the search after waiting for a time ${t}_{\mathrm{wait}}$ post-merger, when ${\dot{f}}_{\mathrm{gw}}|$ decreases such that Equation (13) is satisfied. Alternatively, we can choose shorter ${T}_{\mathrm{drift}}$ (i.e., ${T}_{\mathrm{drift}}\leqslant {\dot{f}}_{\mathrm{gw}}^{-1/2}$) and take ${t}_{\mathrm{wait}}=0$ for all waveforms. However, the sensitivity degrades because the frequency resolution ${\rm{\Delta }}f\gt 1$ Hz is relatively coarse for ${T}_{\mathrm{drift}}\lt 1\,{\rm{s}}$.

In a search without prior knowledge of the signal model, we cover the following parameter space $500\,\mathrm{Hz}\leqslant {f}_{\mathrm{gw}0}\leqslant 2\,\mathrm{kHz}$ for ${10}^{2}\,{\rm{s}}\lesssim \tau \lesssim {10}^{4}\,{\rm{s}}$ using seven discrete ${t}_{\mathrm{wait}}$ values in the range $0\leqslant {t}_{\mathrm{wait}}\leqslant 400\,{\rm{s}}$. Monte-Carlo simulations show that the impact on sensitivity from the mismatch in ${t}_{\mathrm{wait}}$ caused by the granularity is negligible.

We summarize here some of the key technical details concerning the search at hand, while the complete derivation of the search method is in Oliver et al. (2019).

Coherence times. The first step is to select the coherent integration time ${T}_{\mathrm{SFT}}$; i.e., the time-baseline of the SFTs. This cannot be arbitrarily large: to avoid the signal shifting more than half a frequency bin, ${T}_{\mathrm{SFT}}$ must satisfy $| {\dot{f}}_{\mathrm{gw}}| {T}_{\mathrm{SFT}}\,\leqslant 1/(2{T}_{\mathrm{SFT}})$. The time variation of ${f}_{\mathrm{gw}}(t)$ is due to two effects: the spin-down of the source, and the Doppler modulation due to the Earth's motion. It is important to notice that in contrast to the continuous wave case, this method assumes that the Doppler modulation is a subdominant effect. Thus, ${T}_{\mathrm{SFT}}$ can be estimated as:

$$\begin{eqnarray}&&{T}_{\mathrm{SFT}}\leqslant \displaystyle \frac{\sqrt{(n-1)\tau }}{\sqrt{2{f}_{\mathrm{gw}0}}}.\end{eqnarray} \tag{ 18 }$$

Hough transform. Second, each of these SFTs is digitized by setting a threshold ${\rho }_{\mathrm{th}}$ on the normalized power, which is directly related to the false-alarm rate α and false dismissal rate β of the search; the optimal value is 1.6 as derived in Krishnan et al. (2004). This digitized spectrum is then weighted based on the noise floor of the detector and the amplitude modulation of the source. The derivation of the weights is given in Oliver et al. (2019) and in Sintes & Krishnan (2007) for the CW all-sky case.

Detection statistic and significance threshold. Finally, each template—defined by the set of signal parameters $({f}_{\mathrm{gw}0},n,\tau )$—is incoherently integrated through the appropriate summation, known as the number count, over the weighted digitized spectrum following Equation (3). The critical ratio Ψ is defined to evaluate the significance of a given template, based on the results obtained for the weighted number count and its estimates over Gaussian noise for the mean μ and the standard deviation σ:

$$\begin{eqnarray}&&{\rm{\Psi }}=\displaystyle \frac{n-\mu }{\sigma }=\displaystyle \frac{{\displaystyle \sum }_{i=1}^{{N}_{\mathrm{SFTs}}}{w}_{i}{y}_{i}-{\displaystyle \sum }_{i=1}^{{N}_{\mathrm{SFTs}}}{w}_{i}\alpha }{\sqrt{{\displaystyle \sum }_{i=1}^{{N}_{\mathrm{SFTs}}}{w}_{i}^{2}\,\alpha (1-\alpha )}}.\end{eqnarray} \tag{ 19 }$$

Here, y_i corresponds to the ith digitalized bin in a given templated track, N_SFTs is the number of SFTs and the weights are ${w}_{i}\propto {({f}_{\mathrm{gw},i})}^{2m}({a}_{i}^{2}+{b}_{i}^{2})/{S}_{{\rm{n}},i}$, where a_i and b_i are amplitude functions of the antenna pattern found in  (Jaranowski et al. 1998) at the ith time step and ${S}_{{\rm{n}},i}$ is the power spectral density at that given bin. As mentioned in Sintes & Krishnan (2007) and Oliver et al. (2019), any change in the normalization of the weights w_i will not change the resulting sensitivity and will also leave significances, or critical ratios Ψ, unchanged. As shown in Oliver et al. (2019), the critical ratio for a multidetector search can be written as

$$\begin{eqnarray}&&{{\rm{\Psi }}}_{{\rm{m}}}=\displaystyle \sum _{k=1}^{{N}_{\det }}{{\rm{\Psi }}}_{k}{r}_{k}\end{eqnarray} \tag{ 20 }$$

where r_k is each detector's relative contribution ratio—proportional to the number of SFTs, the noise power spectral density, the antenna pattern and the signal amplitude pattern—and N_det is the number of detectors. For detection purposes, a threshold is placed on the two-detector ${{\rm{\Psi }}}_{{\rm{m}}}$ corresponding to a 5σ false-alarm probability for the entire template bank. An additional single-detector threshold, used for a consistency veto step, is extrapolated from the 5σ threshold on ${{\rm{\Psi }}}_{{\rm{m}}}$, to make the veto safer under consideration of the actual differences in r_k for this dataset.

Candidate follow-up. To verify that the 51 outliers found in the initial search step were produced by noise in the detector, and exclude the possibility of having any actual astrophysical signals among them, we performed an additional follow-up step. For each template corresponding to one of the outliers, an analogous analysis was performed but with the template ${f}_{\mathrm{gw}}(t)$ evolution starting instead 1 hr after the merger. From Equation (3), templates starting at merger time have vanishing overlap with these +1 hr delayed versions of themselves. For all outliers, we find that the +1 hr critical ratios are compatible with the results found in the original search within 8%. Given the 5σ false-alarm threshold imposed, the critical ratio for these templates can be considered as stationary non-Gaussian noise with a false dismissal probability of less that 10⁻⁴.

Full details of the adaptation of the FrequencyHough algorithm (Astone et al. 2014) to the case of rapid-spindown post-merger signals are given in Miller et al. (2018). We summarize here some technical details relevant to the search presented in this paper.

Time offsets and search durations. The search is run in three configurations using varying ${T}_{\mathrm{SFT}}=2,4,8\,{\rm{s}}$, covering different observing times: ${\rm{\Delta }}t=1$ hr after merger with ${T}_{\mathrm{SFT}}=2\,{\rm{s}}$ for signals lasting 700–7000 s, ${\rm{\Delta }}t\sim 1.5\mbox{--}3$ hr after merger with ${T}_{\mathrm{SFT}}=4\,{\rm{s}}$ for signals lasting 2000–16000 s, and ${\rm{\Delta }}t\sim 2\mbox{--}7$ hr after merger with ${T}_{\mathrm{SFT}}=8\,{\rm{s}}$ for signals lasting 8000–40000 s. The corresponding end times are set separately for each detector to guarantee that the same effective amount of data is covered even in the presence of gaps; the latest timestamps analyzed for either detector are approximately 3, 8 and 22 hr after merger in the three configurations (see Figure 5).

Zoom In Zoom Out Reset image size

Coordinate transformation. The parameters of the power-law spindown model are transformed to new coordinates, such that in the new space, the behavior of the signal is linear. If we write $k^{\prime} ={(2\pi )}^{n-1}k$, we can rewrite Equation (3) as:

$$\begin{eqnarray}&&{f}_{\mathrm{gw}}(t)=\displaystyle \frac{{f}_{\mathrm{start}}}{{\left[1+k^{\prime} (n-1){f}_{\mathrm{start}}^{n-1}(t-{t}_{\mathrm{ref}})\right]}^{\tfrac{1}{n-1}}}\end{eqnarray} \tag{ 21 }$$

where

$$\begin{eqnarray}&&k^{\prime} =\displaystyle \frac{1}{\tau {f}_{\mathrm{start}}^{n-1}(n-1)}.\end{eqnarray} \tag{ 22 }$$

Then we can make the following change of coordinates:

$$\begin{eqnarray}&&x=\displaystyle \frac{1}{{f}_{\mathrm{gw}}^{n-1}};\quad {x}_{0}=\displaystyle \frac{1}{{f}_{\mathrm{start}}^{n-1}}.\end{eqnarray} \tag{ 23 }$$

Equation (21) becomes the equation of a line:

$$\begin{eqnarray}&&x={x}_{0}+(n-1)k^{\prime} (t-{t}_{0}).\end{eqnarray} \tag{ 24 }$$

Now, we map peaks in the $(t-{t}_{0},{f}_{\mathrm{gw}})$ plane of the detector to lines in the (x₀, k) plane of the source, where both variables in this space relate to the physical parameters of the source.

Grid setup. Our method can be used to determine if a signal is present in the data and estimate its ${f}_{\mathrm{start}}$, ${\dot{f}}_{\mathrm{start}}$ and n. We search across different braking indices with a fine, nonuniform grid determined by ensuring that by stepping from n to n + dn, a signal is confined to one frequency bin. For each braking index, we do a Hough transform and record the most significant candidates in the map. The grid in x₀ is determined by taking the derivative of Equation (23); the grid in k is created by considering a transformation ${f}_{\mathrm{gw}}\to \,{f}_{\mathrm{gw}}+{{df}}_{\mathrm{gw}}=1/{T}_{\mathrm{SFT}}$ and $k\to \,k+{dk}$, then solving for dk imposing that the spindown remains constant.

Both grids are nonuniform and depend on the frequency band and spindown range that we use; however, we over-resolve the grid in x₀ for most frequency bands using the maximum frequency we are analyzing because it is computationally faster and does not result in a sensitivity loss.

The transformation has the disadvantage that it creates nonuniform noise (so peaks at higher frequencies contribute more to the Hough map). We account for this by extending the frequency band that we wish to analyze, and we then only select candidates from the original band.

Coincidence step. Candidates are considered in coincidence between detectors if the Euclidean distance between their recovered parameters x₀ and k is less than 3 bins.

Candidate follow-up. Our candidate follow-up procedure is as follows: we correct for the phase evolution of the candidate recovered in each detector individually. Ideally, if we correct for exactly the right parameters, then we expect a monochromatic signal in the time/frequency peakmap. If we are a bin or so off, then there is a residual spindown or spinup, but the signal is linear. Therefore, we can apply the original FrequencyHough to search for this signal. After applying the original FrequencyHough to the one surviving candidate from this search and performing coincidences again, we found no coincidence, which indicates that the candidate was false. We use the critical ratio as a way to veto candidates, defined as:

$$\begin{eqnarray}&&\mathrm{CR}=\displaystyle \frac{y-\mu }{\sigma }\end{eqnarray} \tag{ 25 }$$

where y is the number count in a given bin in the Hough map, μ is the average number count of the noise, and σ is the standard deviation of the number counts due to noise. We determine if the critical ratio increases in the follow-up, but find that it does not for our one remaining candidate.

Sensitivity estimation procedure. We then computed the strain sensitivity for the different configurations of our search in the following way: the loudest coincident candidate was selected in each 50 Hz band, for n = 5. Its ${f}_{\mathrm{start}}$ and duration were used to inject signals with initial frequencies ranging from ${f}_{\mathrm{start}}$ to ${f}_{\mathrm{start}}+50\,\mathrm{Hz}$, with the highest possible spindown in our configuration ${\dot{f}}_{\mathrm{start}}=1/{T}_{\mathrm{SFT}}^{2}$. Based on our theoretical estimates for sensitivity, the highest initial spindown corresponds to the most conservative result; i.e., the worst case sensitivity. A total of 100 injections were done for each set of parameters; i.e., each point plotted in Figures 2, 3 and 9 and each row in Table 5. Recovery of an injection is defined in the same way as a detection in the real search: within a coincidence distance of 3 bins.