The Adaptive Transient Hough method for long-duration gravitational wave transients

TL;DR

This paper introduces a semi-coherent method utilizing a weighted Hough transform to detect long-duration gravitational wave transients from newborn neutron stars, accounting for rapid spin-down and amplitude evolution.

Contribution

It presents a novel adaptive transient Hough method tailored for intermediate-duration gravitational wave signals, with detailed implementation, statistical analysis, and sensitivity estimation.

Findings

Successfully applied to GW170817 post-merger signals

Verified sensitivity with simulated signals

Demonstrated effectiveness for signals with rapid spin-down

Abstract

This paper describes a new semi-coherent method to search for transient gravitational waves of intermediate duration (hours to days). In order to search for newborn isolated neutron stars with their possibly very rapid spin-down, we model the frequency evolution as a power law. The search uses short Fourier transforms from the output of ground-based gravitational wave detectors and applies a weighted Hough transform, also taking into account the signal's amplitude evolution. We present the technical details for implementing the algorithm, its statistical properties, and a sensitivity estimate. A first example application of this method was in the search for GW170817 post-merger signals, and we verify the estimated sensitivity with simulated signals for this case.

Equations

$$x(t)=n(t)+h(t),$$

$$h(t)=F_{+}(\text{n},\psi ,t)h_{+}(t)+F_{\times }(\text{n},\psi ,t)h_{\times }(t),$$

$$F_{+}(\text{n},\psi ,t)$$

$$F_{\times }(\text{n},\psi ,t)$$

$$h_{+}(t)$$

$$h_{\times }(t)$$

$$A_{+}(t)$$

$$A_{\times }(t)$$

$$h_0(t)=\frac{4{\pi }^{2}G}{c^4}\frac{I_{zz}ϵ}{d}{f}_{\mathrm{gw}}^2(t),$$

$$h_0(t)=\sqrt{\frac{3}{80\pi }}\frac{G}{c^5}\frac{1}{d}\alpha M{R}^3\tilde{J}{f}_{\mathrm{r}}^3(t),$$

$$h_0(t)=A_{\mathrm{m}}{f}_{\mathrm{gw}}^m(t),$$

$$\dot{\Omega }=-\kappa {\Omega }^n,$$

$$\tau =-\frac{{\Omega }_0^{1-n}}{\kappa (1-n)}.$$

$$\hat{f}_{\mathrm{gw}}(t)=\left\{\begin{array}[]{ll}f_{\mathrm{gw},0}\left(\frac{t-T_{0}}{\tau}+1\right)^{\frac{1}{1-n}}&\mathrm{if}\;t\geq T_{0}\,,\\ 0&\mathrm{if}\;t<T_{0}\,,\end{array}\right.$$

$$f_{\mathrm{gw}}(t)=\hat{f}_{\mathrm{gw}}(t)\Big{(}1+\frac{\textbf{v(t)}\cdot\textbf{n}}{c}\Big{)}\,,$$

$$T_{\mathrm{coh}}\le \frac{\sqrt{(n-1)\tau }}{\sqrt{2f_{\mathrm{gw},0}}}.$$

$${\rho }_{i,k}\approx \frac{2|{\tilde{x}}_i[f_k]{|}^2}{T_{\mathrm{SFT}}{S}_{\mathrm{n}}[f_k]},$$

$$\delta n$$

$$\delta \tau$$

$$\delta n$$

$$\delta \tau$$

$${\lambda }_i=\frac{4|\widetilde{h_i}[f_k]{|}^2}{T_{\mathrm{SFT}}{S}_n[f_k]},$$

$$p({\rho }_i|{\lambda }_i)=2{\chi }^2(2{\rho }_i|2,{\lambda }_i)=\exp (-{\rho }_i-\frac{{\lambda }_i}{2})I_0(\sqrt{2{\lambda }_i{\rho }_i}),$$

$$\text{E}[{\rho }_i]$$

$${\sigma }^2[{\rho }_i]$$

$$\alpha ({\rho }_{\mathrm{th}})$$

$${\beta }_{\mathrm{i}}({\rho }_{\mathrm{th}})$$

$$\eta_{i}(\rho_{\mathrm{th}}|\lambda_{i})=\int_{\rho_{\mathrm{th}}}^{\infty}p(\rho|\lambda)d\rho=\alpha\Big{(}1+\frac{\rho_{\mathrm{th}}}{2}\lambda_{i}+O(\lambda_{i}^{2})\Big{)}.$$

$$\nu =\sum _{i=1}^N{w}_i{\nu }_i,$$

$$A=\sum _{i=1}^N{w}_i,$$

$$||w|{|}^2=\sum _{i=1}^N{w}_i^2,$$

$$⟨\nu ⟩$$

$${\sigma }_{\nu }^2$$

$$p(\nu |{\rho }_{\mathrm{th}},\lambda )=\frac{1}{\sqrt{2\pi {\sigma }^2}}{e}^{-(\nu -A\alpha {)}^2/2{\sigma }^2}.$$

$$\alpha_{\mathrm{I}}=\int_{\nu_{\mathrm{th}}}^{\infty}p(\nu|\rho_{\mathrm{th}},0)d\nu=\frac{1}{2}\mathrm{erfc}\Big{(}\frac{\nu_{\mathrm{th}}-\langle\nu\rangle}{\sqrt{2}\sigma_{\mathrm{\nu}}}\Big{)}\,.$$

Full text

I INTRODUCTION

The advanced gravitational wave (GW) detector era has provided us with multiple detections from binary compact objects (Abbott et al., 2018a) including GW170817, the first binary neutron star (BNS) coalescence (Abbott et al., 2017). This detection motivated the development of the new search method presented in this paper, focusing on the possible birth of a rapidly rotating highly magnetized neutron star (NS) spinning down through some combination of GW and electromagnetic emission. For a very massive remnant, the collapse would occur in a short time scale (as explored in Abbott et al. (2017, 2018b)), but for low total mass and some equations of state, the emitted GW signal could have an intermediate duration on the order of hours to days (Baiotti and Rezzolla, 2017; Piro et al., 2017).

This regime of GW signal durations has long been mostly unexplored from the data analysis side. The expected rapid frequency and amplitude evolution, in combination with observation times still much longer than e.g. for individual binary coalescences, pose unique challenges on analysis algorithms. Other pre-existing or recently developed methods to search for intermediate-duration signals include the Stochastic Transient Analysis Multi-detector Pipeline (STAMP) (Thrane et al., 2011), the Hidden Markov model Viterbi algorithm (Sun and Melatos, 2018) and a generalization of the FrequencyHough method (Miller et al., 2018). The first two are generic unmodeled searches, while the last is a modeled search for power-law spin-downs similar to the one described in this paper. Together with those three pipelines, our new Adaptive Transient Hough (ATrHough) method has already contributed to the search for a long-duration transient signal from a putative NS remnant of GW170817 described in Abbott et al. (2018c).

The Adaptive Transient Hough is a semi-coherent analysis adapted from the SkyHough (Sintes and Krishnan, 2007, 2006; Krishnan et al., 2004) search for continuous wave (CW) signals. Like most other CW searches (Riles, 2017), the original SkyHough assumes a constant intrinsic amplitude and slowly evolving frequency, and hence cannot be used to search for transient GWs with rapid frequency and amplitude evolution (see quantitative comparison in Sec. II), for which we have now specifically developed the new method.

The ATrHough method will also have wider applicability beyond the case of BNS remnants, as signals with similar durations and evolutionary behaviour are also possible from young NSs born through the regular supernova channel (Palomba, 2001; Dall’Osso et al., 2009, 2015; Lasky and Glampedakis, 2016; Dall’Osso et al., 2018), emitted either by r-mode oscillations (Owen et al., 1998; Andersson and Kokkotas, 2001) or quadrupolar deformations.

The paper is organized as follows: section II briefly describes the expected signal from a remnant NS. Section III summarizes the general strategy of a hierarchical search and its implementation, section IV studies its statistical properties, and section V introduces the threshold and vetoes required for a robust detection strategy. Finally section VI presents an estimate for the search sensitivity and section VII presents our conclusions.

II THE TRANSIENT SIGNAL MODEL

The output of a GW detector can be represented by

[TABLE]

where $n(t)$ is the detector noise at time $t$, and $h(t)$ is the strain induced by a GW signal:

[TABLE]

where $F_{+,\times}$ are the detector antenna patterns, which depend on a unit-vector n corresponding to the sky location of the source and on the wave polarization angle $\psi$, and vary with time due to the movement of the detector frames with the Earth. For ground-based detectors with perpendicular arms, the expressions for $F_{+,\times}$ are (Jaranowski et al., 1998):

[TABLE]

where the functions $a(t;\textbf{n})$ and $b(t;\textbf{n})$ are independent of $\psi$. For convenience, we do not explicitly write out the n and $\psi$ dependence from here on. Now the waveforms for the two polarizations $h_{+,\times}$ are:

[TABLE]

where $\Phi(t)$ is the phase evolution of the signal and $A_{+,\times}(t)$ are the (time-varying) amplitude parameters depending on the orientation $\cos{\iota}$ of the source and on the strain amplitude evolution $h_{0}(t)$ as follows:

[TABLE]

The time evolution of the dimensionless strain amplitude $h_{0}(t)$ depends on the emission mechanism; if for example it is due to a constant non-axisymmetrical deformation of the source NS, but the frequency decays over time, the amplitude evolves as

[TABLE]

where $c$ is the speed of light, $I_{zz}$ is the z-z component of the star’s moment of inertia with the z-axis being its spin axis, $\epsilon:=(I_{xx}-I_{yy})/I_{zz}$ is the equatorial ellipticity of the star, and $d$ is its distance from Earth. Another mechanism covered by this method is GW emission from r-mode oscillations, which are the result of small velocity and density perturbations of the NS fluid, causing a time-varying moment of inertia restored throw Coriolis force; for these, the amplitude evolution is given by

[TABLE]

where $f_{\mathrm{r}}$ is the rotation frequency of the source, $\widetilde{J}=0.01635$ is a dimensionless constant, $M$ is the NS mass, $R$ its radius and $\alpha$ is a dimensionless amplitude described in more detail in Owen et al. (1998).

Independent of the specific emission scenario, the amplitude evolution $h_{0}(t)$ can be written in a more general form as:

[TABLE]

where $m$ and $A_{\mathrm{m}}$ are constants defined by the emission mechanism.

To characterize the frequency evolution of a newborn NS we apply the waveform model from Lasky et al. (2017); Sarin et al. (2018), originating from the general torque equation

[TABLE]

where $\Omega$ and $\dot{\Omega}$ are the frequency of rotation of the source and its derivative. (When we focus on GW emission due to a non-axisymmetrical shape and do not consider the free precession case (Zimmermann and Szedenits, 1979; Jones and Andersson, 2001), the frequency of GW emission is $f_{\mathrm{gw}}=\Omega/\pi$.) Furthermore, $n$ is called the star’s braking index and $\kappa$ is associated to the spindown timescale:

[TABLE]

The solution of Eq. (9) for arbitrary braking index $n$ characterizes the frequency evolution:

[TABLE]

where $f_{\mathrm{gw},0}$ corresponds to the frequency at the start of the emission ($t=T_{0}$); for simplicity let us set $T_{0}=0$ s. A braking index of $n=5$ corresponds to pure GW emission from a non-axisymmetric rotator. This equation can also be applied to r-modes, for which $n\lesssim 7$ (Alford and Schwenzer, 2014, 2015).

The Eq. (11) frequency evolution model and resulting amplitude evolution as per Eq. (8) is the key difference between our new search method and the SkyHough search (Krishnan et al., 2004) for CW signals, which instead uses a Taylor expansion for the slowly-evolving frequency of mature NSs and assumes constant intrinsic amplitude.

To demonstrate explicitly that such an expansion is unsuited to search for signals with rapid spindowns, let us consider that the frequency resolution of a fully-coherent CW-like search over an observation time is $\delta f_{\mathrm{gw}}=1/T_{\mathrm{obs}}$. Hence, for a Taylor expansion model $\mathcal{T}[f_{\mathrm{gw}}(t),s]$ to order $s$ in $f_{\mathrm{gw}}(t)$, the requirement is $|f_{\mathrm{gw}}(t)-\mathcal{T}[f_{\mathrm{gw}}(t),s]|<1/T_{\mathrm{obs}}$. Now we see that at least a 16th order expansion is required to track sources with astrophysically relevant example parameters (compare Abbott et al. (2018c)) $f_{\mathrm{gw}}(0)=1000$ Hz, $\tau=10000$ s and $n=5$ over $T_{\mathrm{obs}}=5000$ s, making this approach computationally prohibitive. On the other hand, the search method introduced in the following uses the exact analytical form with its only three free parameters $(n,f_{\mathrm{gw},0},\tau)$ to create a template grid that ensures complete coverage, while keeping the analysis computationally feasible.

As in other semi-coherent searches, this method considers as negligible – and therefore ignores – relativistic corrections, and those due to the time delay between the detector and the solar-system barycenter (SSB). Therefore only the instantaneous signal frequency in the detector frame needs to be calculated:

[TABLE]

where v(t) is the detector velocity with respect to the SSB frame. Note that now the time coordinate $t$ corresponds to time at the detector.

III The Adaptive Transient Hough Method

This section discusses the implementation of the Adaptative Transient Hough (ATrHough) method, a pipeline based on the semi-coherent SkyHough search for CWs described in Krishnan et al. (2004); Sintes and Krishnan (2007). The common ground of both searches is the use of a weighted Hough transform on Short-time Fourier Transforms (SFTs) as the input data. The Hough transform is an algorithm widely used in pattern recognition; here the pattern is defined by the frequency evolution of the signal in the detector data. In both CW and transient cases, the weights take into account the amplitude modulation of the signal, caused by the antenna pattern, and the changing noise floor between SFTs. But as a difference to the CW SkyHough search, the new ATrHough method also includes the source amplitude evolution in the weights.

Together with the power-law frequency evolution model from Eq. (11), the amplitude weights allow a sensitive search for transient signals from rapidly evolving newborn NSs. Meanwhile, the main framework and statistical properties are the same as in the SkyHough method. In the following we summarize them in the new context, and add the required transient-specific details.

III.1 Length of Short-duration Fourier Transforms

We first obtain a collection of SFTs by dividing the full observation time $T_{\mathrm{obs}}$ in $N$ segments of length $T_{\mathrm{coh}}$. The maximum length of $T_{\mathrm{coh}}$ is calculated by imposing the 1/4-cycle criterion introduced in Jaranowski et al. (1998): This leads to a requirement $2|df/dt|\leq T_{\mathrm{coh}}^{-2}$, ensuring that the maximum modulation corresponds to only half a bin at the search frequency resolution $\delta f=1/T_{\mathrm{obs}}$. From Eq. (12) the spin-down modulation is given by two effects, the spin-down of the source and the Doppler modulation resulting from the Earth’s motion. The constraint imposed by the spin-down of the source is:

[TABLE]

The range of maximum allowed $T_{\mathrm{coh}}$ for the parameter space covered in Abbott et al. (2018c) is on the order of seconds, as shown in Fig. 1. On the other hand, the constraint imposed by Doppler modulation is on the order of hours, as discussed in Krishnan et al. (2004). Therefore we will consider only the spin-down of the source as the dominant threshold for $T_{\mathrm{coh}}$.

III.2 The peak-gram

The Hough transform requires a digitized spectrum as its input, with time-frequency bins categorized in two classes. The ATrHough generates this by setting a threshold $\rho_{\mathrm{th}}$ on the normalized power spectrum $\rho_{i}$ to conduct the bin selection:

[TABLE]

where $[.]$ indicates a discrete series and the index $i$ corresponds to the $i^{\mathrm{th}}$ time step. That is, $\tilde{x}_{i}[f_{k}]$ is the value obtained from the $i^{\mathrm{th}}$ SFT on the $k^{\mathrm{th}}$ frequency bin. Furthermore, $S_{\mathrm{n}}$ is the single-sided Power Spectral Density (PSD) of the noise in the same bin. In the following, we drop the explicit $k$ index, as we are only interested in the frequency bins following the signal track. If $\rho_{i}\geq\rho_{\mathrm{th}}$, then a value of $1$ is assigned to that bin, and a [math] otherwise. The result of this process is known as the peak-gram.

III.3 Resolution in $\tau$ and $n$ space

The Hough transform is applied to find the statistical significance of each template in a bank over parameter space. A template is defined by the intrinsic parameters of the signal, $\vec{\xi}=(f_{\mathrm{gw},0},n,\tau,T_{0})$. To conduct a wide-parameter space search, we create a grid that ensures contiguous templates to deviate from each other by at most one frequency bin over a duration $T_{\mathrm{obs}}$; this ensures the computation of at least all independent templates (by the 1/4-cycle criterion) between $t=0$ s and $t=T_{\mathrm{obs}}$. The grid is constructed with the following step sizes:

[TABLE]

where $\delta f=1/T_{\mathrm{coh}}$. Hence,

[TABLE]

The two grid step sizes are inversely proportional to $f_{\mathrm{gw},0}$. Fig. 2 represents the obtained $\delta\tau$ and $\delta n$ for a fixed $T_{\mathrm{coh}}$, $T_{\mathrm{obs}}$ and $f_{\mathrm{gw},0}$ inside the $\tau$, $n$ ranges.

The practical implementation of the grid is defined by a nested loop; a pipeline diagram can be seen in Fig. 3. First, we select the minimum value of $\delta n$ over the $\tau$ range as shown in Fig. 4, given a set of $(T_{\mathrm{obs}},T_{\mathrm{coh}},n)$ and the maximum $f_{\mathrm{gw},0}$; then we calculate $\delta\tau$ as in Fig. 5. We will recalculate $\delta n$ and $\delta\tau$ on each iteration of the $n$ and $\tau$ loops respectively.

In order to reduce the number of templates or grid points required by the search, we need to split the $\tau$ and $f_{\mathrm{gw},0}$ ranges of the whole search space into smaller subdomains. To do so, we will typically create bands for $\tau$ smaller than $10\%$ of $T_{\mathrm{obs}}$ and frequency bands between 50 and 100 Hz in width. Each sub-domain will be analyzed independently, making the computational load smaller. It is possible to make the domains larger, but the necessary refinement of the grid in certain areas will make the search less computationally efficient overall.

Fig. 6 shows the distribution and number of templates used for different $T_{\mathrm{obs}}$ given a search that covers an analogous parameter space as Abbott et al. (2018c). Here templates are calculated with the maximum integer coherence length allowed, and the minimum $T_{\mathrm{coh}}$ considered for this figure and the search is 1 s.

IV STATISTICAL PROPERTIES

IV.1 The coherent statistic

For the following section we make the assumption of stationary Gaussian noise with zero mean in order to characterize the output of the detectors, for which the normalized power $2\rho_{i}$ in the presence of a signal $h$ follows a non-central $\chi^{2}$ distribution with 2 degrees of freedom and a non-centrality parameter

[TABLE]

where $|\widetilde{h_{i}}[f_{k}]|$ is the Fourier transform of the signal and, as before in Eq. (14) for the normalized power $\rho_{i}$, for $\lambda_{i}$ we suppress the $k$ dependence. Then the probability distribution for $\rho_{i}$ is:

[TABLE]

where $I_{0}$ is the zero-order modified Bessel function of the first kind.

The mean and variance for this distribution are respectively:

[TABLE]

The false alarm and false dismissal probabilities for a frequency bin to be above the power spectrum threshold are:

[TABLE]

The probability $\eta_{i}$ that a given frequency bin is selected is, in the small-signal approximation:

[TABLE]

IV.2 The incoherent number-count statistic

If a signal is present, the non-centrality parameter $\lambda_{i}$ will change for different SFTs. As pointed out previously, this can happen both because the noise may not be stationary and because the amplitude modulation of the signal changes over time. In other words, the observed signal power $|h|^{2}$ changes due to the non-uniform antenna pattern of the detector and due to the intrinsic spindown. Therefore, the detection probability $\eta_{i}$ changes across SFTs. This is taken into account by adapting the non-demodulated weighted Hough approach mentioned before and covered in Sintes and Krishnan (2007); it is a similar strategy to the one applied in the StackSlide (Brady and Creighton, 2000) and PowerFlux (Dergachev, 2005; Dergachev and Riles, 2005) algorithms. The starting point is to generalize the integer number-count statistic, which we would obtain directly from the peak-map, to a non-integer weighted statistic

[TABLE]

where $N$ is the number of SFTs, $\nu_{i}$ is the value assigned to the bin selected from the peak-gram in the $i^{\mathrm{th}}$ time step for the current template, and $w_{i}$ are a constant set of weights given for each template with $w_{I}\propto 1/N$. It is important to notice that in order to maximize the sensitivity of the search the selection of weights is not arbitrary; we will derive the optimal choice in Sec. IV.4. For now, we define the normalization terms

[TABLE]

This step in the search (computing $\nu$) is known as the incoherent sum; the templates in a search are then ranked based on their number count $\nu$. Applying the linearity of the expectation value, the mean and variance for the incoherent step in the absence of a signal are:

[TABLE]

As shown in Sintes and Krishnan (2007) and applied in multiple CW searches like Astone et al. (2014), when optimal weights are chosen we can, for a sufficient number of SFTs, evaluate the significance of an observation by approximating the number count distribution by a Gaussian with the right mean and variance:

[TABLE]

This becomes a very good approximation for $N>1000$, and e.g. the typical number of SFTs searched in Abbott et al. (2018c) is indeed above that number. We provide some empirical tests of this approximation in appendix A.

Thus one can derive the number count threshold $\nu_{th}$ based on the incoherent false-alarm rate as

[TABLE]

For a given set of weights and peak selection threshold, this equation decides what number count threshold must be used to obtain a desired $\alpha_{\mathrm{I}}$. We can solve this as

[TABLE]

The false-dismissal rate requires the computation of the mean and variance, which in the presence of a small signal are:

[TABLE]

If the small-signal approximation is applied, $\sigma_{\mathrm{\eta}}^{2}$ can be expanded to first order in $\lambda_{i}$:

[TABLE]

We again approximate the number count distribution $p(\eta|h)$ by a Gaussian distribution with the above mean and variance, yielding the false-dismissal rate as follows:

[TABLE]

IV.3 Setting up the threshold

Considering the statistical significance in a template as $\,s:=1-\alpha_{\mathrm{I}}-\beta_{\mathrm{I}}$ and using the properties of the complementary error function, we can introduce a quantity

[TABLE]

This equation can be shown to reduce to $s$ when $S=0$, and as it grows monotonically we can take it as a measure of the statistical significance of the search. By expanding to the first order in $\lambda_{i}$, we derive the following expression:

[TABLE]

Imposing again optimal weights which are proportional to $1/N$, for large values of $N$ the first term on the right-hand side of this equation is proportional to $\sqrt{N}$, while the second term does not grow with $N$. Thus the first term dominates, yielding

[TABLE]

The peak selection threshold is chosen to minimize $\beta_{\mathrm{I}}$, or equivalently maximize $S$ for fixed $\alpha_{\mathrm{I}}$:

[TABLE]

As derived in Krishnan et al. (2004), this threshold is independent of the choice of weights; and the solution of the previous equation is $\rho_{\mathrm{th}}=1.6$ which leads to $\alpha=e^{-\rho_{\mathrm{th}}}=0.2$. Different thresholds can be imposed, yielding different $\alpha$, but they would not maximize the statistical significance of the template.

IV.4 Calibration of the weights

To define an appropriate set of weights, we start by considering the modulus square of the signal’s Fourier transform on the $i^{\mathrm{th}}$ SFT, depending on the antenna patterns $F_{+,\times}$ from Eq. (3) and amplitudes $A_{+,\times}$ from Eq. (5):

[TABLE]

From here on, the subindex $i$ runs over segments and in the case of a fuction it imposes a time average of length $T_{\mathrm{coh}}$, e.g. for the time-evolving GW frequency from Eq. (11): $f_{\mathrm{gw},i}=\int_{T_{i}-T_{\mathrm{coh}}/2}^{T_{i}+T_{\mathrm{coh}}/2}f_{\mathrm{gw}}(t)dt/T_{\mathrm{coh}}$. The subindex $k$ corresponds to the $k^{\mathrm{th}}$ frequency bin, selected so that $f_{\mathrm{gw},i}\in(f_{k}-\delta f/2,f_{k}+\delta f/2)$. The average over that interval is

[TABLE]

Now we can average over the NS’s orientation $\cos\iota$ and the polarization angle $\psi$ appearing in the antenna patterns and find the following relationships:

[TABLE]

where $h_{0,0}=h_{0}(t=t_{0})$ is the initial amplitude at $t_{0}$.

Combining all these results:

[TABLE]

and substituting this into Eq. (33), the sensitivity is

[TABLE]

From this, we see that the sensitivity is related to the detector response and the amplitude modulation of the signal, which we can summarize in a quantity

[TABLE]

Calculating the maximum of the inner product $\textbf{w}\cdot\textbf{X}$ shows that the weights guarantee the best sensitivity for a given template if the two vectors are proportional to each other, i.e. $w_{i}\propto X_{i}$. At the same time, we see that any overall rescaling of the weights ($\hat{w}_{i}=kw_{i}$) has no impact on $S$, as for any constant $k$ the value of the detectable dimensionless strain amplitude $h_{\mathrm{0,0}}$ at $t=0$ s remains unchanged.

In summary, as also illustrated for an example simulated signal in Fig. 7, the use of appropriate weights ensures our search properly accounts for both the source’s amplitude decay and the effects of the detector response changing with time. This gives us the ability to compare templates across the search parameter space, comparing very fast frequency decays with slower ones.

If the value $\rho_{\mathrm{th}}$ = 1.6 is substituted in Eq. (39), the minimum theoretical value of $h_{0}$ that the search can recover is:

[TABLE]

IV.5 Critical Ratio $\Psi$

The critical ratio $\Psi$ is a new statistic that quantifies the significance of a given template. Based on the weighted number count and quantities from Eqs. (22)–(24), we define

[TABLE]

As mentioned before, any normalization of the weights will not change the sensitivity of the search. It will also leave the significance or critical ratio in each template unchanged. Considering the previous equation as the single-detector case, the multi-detector critical ratio is defined as

[TABLE]

where $N_{\mathrm{M}}$ is the number of detectors and $N_{k}$ is the number of SFTs in detector $k$, while $w_{i;k}$ and $\nu_{i;k}$ are the weights and number count assigned to the $i^{th}$ SFT for that detector and a given template. We can also rewrite this as

[TABLE]

where $\Psi_{k}$ is the critical ratio for each single detector $k$.

In a multi-detector search, the duty factors (fraction of time a detector is recording usable data) and noise floors may differ between detectors. To quantify the contribution of each detector to the multi-detector critical ratio, the relative contribution ratio is defined as

[TABLE]

Using the previous equations, the critical ratio for a multi-detector search takes a very simple form:

[TABLE]

V Vetoes on Critical ratio and Time Consistency

Candidates that appear significant by their critical ratio can be due to astrophysical sources, but also due to non-Gaussian noise artifacts in the data. To make the search robust against such artifacts, we introduce vetoes that test for each candidate (i) its consistency between detectors and (ii) the consistency of its transient behavior with the target astrophysical model.

V.1 The Critical ratio $\Psi$-veto

The threshold for a search is determined under the assumption of detector noise following a stationary zero-mean Gaussian distribution with a power spectral density $S_{\mathrm{n}}(f)$. A template is considered as a candidate when its $\Psi$ exceeds a pre-specified threshold for which the probability of a false alarm due to noise alone is small. (See Fig. 8.) The overall false-alarm probability $\alpha_{\mathrm{S}}$ of the search can be approximated as the product of the number of trials (i.e number of templates $N_{\mathrm{t}}$) and the previously introduced false-alarm probability $\alpha_{\mathrm{I}}$. Now we can rewrite Eq. (26) in terms of the critical-ratio threshold $\Psi_{\mathrm{th}}$:

[TABLE]

If the critical ratio in a template exceeds the threshold, as a follow-up veto we can rephrase the question and consider each detector as an independent single trial, obtaining a threshold $\Psi_{\mathrm{th}}^{\mathrm{D}}$ for each detector. This threshold will correspond to Eq. (47) with $N_{\mathrm{t}}=1$ and any given template that fails to satisfy it in either detector will be vetoed.

V.2 The time-inconsistency veto

To check that the transient behavior of the signal matches our model, we introduce an additional veto. Let us consider a candidate template $\vec{\xi_{\mathrm{C}}}=(f_{\mathrm{gw},0},n,\tau,T_{0}=T_{\mathrm{event}})$ and a time-shifted version $\vec{\xi_{\mathrm{F}}}=(f_{\mathrm{gw},0},n,\tau,T_{0}=T_{\mathrm{event}}+T_{\mathrm{F}})$. These will be completely independent if $T_{\mathrm{F}}=-T_{\mathrm{obs}}$; see Fig. 9 for an example. Other time shifts could be used for a veto as well, as long as the contribution of the candidate signal $\vec{\xi_{\mathrm{C}}}$ to $\Psi_{\mathrm{F}}$ of the shifted template $\vec{\xi_{\mathrm{F}}}$ is zero.

The obtained value $\Psi_{\mathrm{F}}$ will indicate how much of the original candidate’s $\Psi_{\mathrm{C}}$ seems to come from a stationary contribution instead. Stationary spectral line artifacts are common in the LIGO data (Covas et al., 2018) and hence this veto is important to remove non-astrophysical false candidates. In other words, we assign a probability to stationary lines to be the cause of the candidate. To estimate this probability we reuse Eq. (47) for a single follow-up trial. If the resulting probability corresponds to more than 6 sigmas, we can safely reject the candidate.

VI Search sensitivity

In Eq. (41) we have obtained an estimate for the sensitivity of a search as the smallest amplitude that would cross the number-count threshold for a given false-alarm rate $\alpha_{\mathrm{I}}$ and false-dismissal rate $\beta_{\mathrm{I}}$.

As a specific astrophysical case, let us concentrate on the isolated non-axisymmetric magnetar scenario as considered in the GW170817 long-duration postmerger search (Abbott et al., 2018c). In this model, the amplitude exponent $m$ in Eq. (8) takes a nominal value of 2 and the signal amplitude $h_{0}(t)$ is given by Eq. (6).111In the case of GWs emitted from r-mode oscillations, we have instead $n\lesssim 7$, $m=3$ and $h_{0}(t)$ is given by Eq. (7). This case is described in more detail e.g. in Owen et al. (1998).

In Fig. 10 we show an example signal recovery for the same injection as in Fig. 7. As we can see, power-law templates in principle allow to succesfully track this type of signal even without weights, but including the source’s amplitude decay in the weights from Sec. IV.4 is crucial for robust recovery and to fully profit from long observation times.

Combining the amplitude from Eq. (6) with the sensitivity as given by Eq. (41), the astrophysical range of the search is

[TABLE]

We now calculate an astrophysical range estimate for a search setup corresponding to the ATrHough analysis performed as one of four searches in Abbott et al. (2018c). We use the aLIGO O2 sensitivity $S_{\mathrm{n}}$ during the GW170817 event to calculate the weights, and for the remnant’s moment of inertia we use the same value as in Abbott et al. (2018c), $I_{\mathrm{zz}}=100M_{\odot}^{3}G^{2}/c^{4}\approx 4.34\times 10^{38}\,\mathrm{kg}\,\mathrm{m}^{2}$.

In Fig. 11 we compare the analytical estimate with the empirical recovery fraction for a set of injections. Those were originally performed for the sensitivity estimate in the GW170817 post-merger search (Abbott et al., 2018c). The recovery criterion corresponds to $\Psi_{\mathrm{th}}=9$, or a $5\sigma$ significance. We have concentrated here on a braking index $n=5$ that corresponds to pure GW emission, and covered ranges of $f_{\mathrm{gw},0}$ and $\tau$ as shown in Fig. 12. The procedure to obtain the experimental results consisted in selecting 10 Hz wide frequency bands, for each band injecting 1000 simulated signals into O2 data with amplitudes around the astrophysical range estimate. The purpose was to find the amplitude corresponding to $90\%$ recovery efficiency. The parameters $\tau$ and $f_{\mathrm{gw},0}$ were randomized within 10 bins of their nominal value; i.e. the injection parameters are not perfectly aligned with the search grid, thus allowing for a realistic exploration of search mismatch in the recovery.

We do not expect an exact agreement between analytical prediction and sensitivity measured from injections, as the analytical estimate is based on a Gaussian noise approximation. But the results are sufficiently close to demonstrate that Eq. (41) is useful for the purpose of setting up future searches.

VII CONCLUSIONS

In this paper we have described a new semi-coherent search method for quasi-monochromatic gravitational waves, using short incoherent steps of the order of seconds with the intention to track signals of intermediate durations (of the order of hours to days) even if these show rapid frequency evolution. The main innovations compared to previous versions of the Hough transform method (Sintes and Krishnan, 2007, 2006; Krishnan et al., 2004) are the new frequency-evolution templates and the additional inclusion of amplitude evolution in the Hough weights.

In introducing this new method and estimating its sensitivity, we have concentrated on the model of power-law spin-down for a newborn NS. As applied in the GW170817 post-merger remnant search (Abbott et al., 2018c), the astrophysical range of this method at 90% detection confidence is at $\sim 1$ Mpc with LIGO sensitivity at the end of the second observing run. With future instruments like the Einstein Telescope (Punturo et al., 2010; Hild et al., 2011; Sathyaprakash et al., 2012), this range could increase by a factor of $\sim 20$.

One disadvantage of modeled semi-coherent methods like this one is the need to explicitly set a starting time for the signal model. On the other hand, it is a suitable method to perform fast and economic follow-ups of known merger events or for promising candidates identified by more generic searches, allowing to reliably set up a fixed false-alarm rate of the overall search.

The same strategy can also easily be translated to signals following other spin-down patterns than the power-law model we focused on so far, with the definition of weights and parameter space grids following the same procedure as introduced in this paper.

Acknowledgments

We thank the LIGO-Virgo Continuous Wave working group and the GW170817 postmerger search team, in particular S. Banagiri, M. Bejger, A. Miller, L. Sun, K. Wette and S. Zhu, for many fruitful discussions. M.O. and A.M.S. acknowledge the support of the Spanish Agencia Estatal de Investigación and Ministerio de Ciencia, Innovación y Universidades grants FPA2016-76821-P, FPA2017-90687-REDC, FPA2017-90566-REDC, and FPA2015-68783-REDT, the Vicepresidencia i Conselleria d’Innovació, Recerca i Turisme del Govern de les Illes Balears, the European Union FEDER funds, and the EU COST actions CA16104, CA16214 and CA17137. The authors are grateful for computational resources provided by the LIGO Laboratory and supported by National Science Foundation Grants PHY-0757058 and PHY-0823459.

Appendix A Testing the Gaussian approximation for the weighted number count $\nu$

In Eq. (25) we have approximated the distribution $p(\nu|\rho_{\mathrm{th}},\lambda)$ of the weighted number-count statistic $\nu$, when using appropriate weights and for a sufficient number of SFTs, by a Gaussian. Here we present some simple empirical tests of this limiting behaviour in configurations similar to the search implemented in Abbott et al. (2018c).

Using the same machinery as before, we have analysed 100 simulated data sets, each consisting of 1000 segments of Gaussian noise with no GW injection ($h_{0}=0$). For each, we have computed the number count for 10000 template trials, covering a small fraction of the parameter space around a random point corresponding to the ‘null injection’, and using the weights proportional to $X_{i}$ as introduced in Sec. IV.4. We have then compared the resulting empirical distribution of $\nu$ with our Gaussian approximation from Eq. (25). An example is shown in Fig. 13 to illustrate the agreement between the two distributions.

To further evaluate the (dis-)agreement between two distributions $P$ and $Q$, one can compute the Kullback-Leibler (KL) divergence (Kullback and Leibler, 1951) (in bits):

[TABLE]

for a discrete set $\mathcal{X}$ of measured values. Note the asymmetry in this definition; here we take the Gaussian for $P$ and the empirical results for $Q$. A histogram of KL divergences between the Gaussian from Eq. (25) and the empirical distributions from the 1000 simulations is shown in Fig. 14. We see that there is far less than 1 bit of information between the two distributions in all draws. Hence, based on the KL divergence we can consider the approximation from Eq. (25) as a sufficiently robust basis for estimating significance of our search results.