Application of a Hough search for continuous gravitational waves on data from the fifth LIGO science run

The focus of this article is the search for evidence of continuous gravitational waves (GWs), as might be radiated by nearby, rapidly spinning neutron stars, in data from the Laser interferometer gravitational-wave observatory (LIGO) [1]. The data used in this paper were produced during LIGO's fifth science run (S5) that started on November 4, 2005 and ended on October 1, 2007.

Spinning neutron stars are promising sources of GW signals in the LIGO frequency band. These objects may generate continuous GWs through a variety of mechanisms including non-axisymmetric distortions of the neutron star, unstable oscillation modes in the fluid part of the star and free precession [2–6]. Independently of the specific mechanism, the emitted signal is a quasi-periodic wave whose frequency changes slowly during the observation time due to energy loss through GW emission, and possibly other mechanisms. At an Earth-based detector the signal exhibits amplitude and phase modulations due to the motion of the Earth with respect to the source.

A number of searches have been carried out previously in LIGO data [7–18] including: targeted searches in which precise pulsar ephemerides from radio, x-ray or γ-ray observations can be used in a coherent integration over the full observation span; directed searches in which the direction of the source is known precisely, but for which little or no frequency information is known; and all-sky searches in which there is no information about location or frequency.

All-sky searches for unknown neutron stars must cope with a very large parameter space volume. Optimal methods based on coherent integration over the full observation time are completely unfeasible since the template bank spacing decreases dramatically with observation time, and even for a coherent time baseline of just few days, a wide-frequency-band all-sky search is computationally extremely challenging. Therefore hierarchical approaches have been proposed [19–24] which incorporate semi-coherent methods into the analysis. These techniques are less sensitive for the same observation time but are computationally inexpensive. The Hough transform [7, 10, 21, 25, 26] is an example of such a method and has been used in previous wide-parameter-space searches published by the LIGO and Virgo Collaborations. Moreover it has also been used in the hierarchical approach for Einstein@Home searches, as the incoherent method to combine the information from coherently analyzed segments [14, 18].

In this paper we report the results of an all-sky search making use of the 'weighted Hough' method [10, 25, 26]. The 'weighted Hough' was developed to improve the sensitivity of the 'standard Hough' search [7, 21] and allows us to analyze data from multiple detectors, taking into account the different sensitivities.

The work presented here achieves improved sensitivity compared to previous Hough searches [7, 10] by splitting the run into two year-long portions and requiring consistency between signal levels in the two separate years for each candidate event, in addition to incorporating a χ²-test [27]. This new pipeline is efficient at rejecting background, allowing us to lower the event threshold and achieve improved sensitivity. The parameter space searched in our analysis covers the frequency range 50 < f < 1000 Hz and the frequency time-derivative range $-8.9\times 10^{-10}<\dot{f}<0$ Hz s⁻¹. We detect no signals, so our results are presented as strain amplitudes h₀ excluded at 95% confidence, marginalized over the above $\dot{f}$ interval.

Through the use of significant distributed computing resources [28], another search [18] has achieved better sensitivity on the same data as the search described here. But the Einstein@Home production run on the second year of S5 LIGO data required about 9.5 months, used a total of approximately 25 000 CPU (central processing unit) years [18], and required five weeks for the post-processing on a cluster with 6720 CPU cores. The search presented in this paper used only 500 CPU months to process each of the two years of data, representing a computational cost more than two orders of magnitude smaller. This is also an order of magnitude smaller than the computational cost of the semi-coherent 'PowerFlux' search reported in a previous paper [17].

The significance of our analysis is through offering an independent analysis to cross-check these results, and a method that allows the attainment of sensitivity close to that of the Einstein@Home search at substantially reduced computational burden. This technique will be particularly important in the advanced LIGO and Virgo detector when applied to 'quick-look' searches for nearby sources that may have detectable electromagnetic counterparts. Moreover, the Hough transform is more robust than other computationally efficient semi-coherent methods with respect to noise spectral disturbances [10] and phase modeling of the signal. In particular, it is also more robust than the Einstein@Home search to the non inclusion of second order frequency derivatives.

An important feature to note is that the sensitivity of the Hough search is proportional to $1/( N^{1/4} \sqrt{T_{{{{\rm coh}}}}})$ or $N^{1/4}/ \sqrt{T_{{{{\rm obs}}}}}$, assuming T_obs = NT_coh, being N the number of data segments coherently integrated over a time baseline T_coh and combined using the Hough transform over the whole observation time T_obs, while for a coherent search over the whole observation time, the sensitivity is proportional to $1/ \sqrt{T_{{{{\rm obs}}}}}$. This illustrates the lost of sensitivity introduced combining the different data segments incoherently but, of course, this is compensated by the lesser computational requirements of the semi-coherent method.

For sufficiently short segments (T_coh of the order of 30 min or less), the signal remains within a single Fourier frequency bin in each segment. In this case a simple Fourier transform can be applied as a coherent integration method. As the segment duration T_coh is increased, it becomes necessary to account for signal modulations within each segment by computing the so-called $\mathcal{F}$-statistic [29] over a grid in the space of phase evolution parameters, whose spacing decreases dramatically with time baseline T_coh. This results in a significant increase in the computational requirements of the search and also limits the significant thresholds for data points selection and the ultimate sensitivity of the search.

The search presented here is based on 30 min long coherent integration times, being this the reason for the significant reduction of the computational time compared to the Einstein@Home search [18] in which the span of each segment was set equal to 25 h. For an in-depth discussion on how to estimate and optimize the sensitivity of wide area searches for spinning neutron stars at a given computational cost, we refer the reader to [22, 23, 30].

This paper is organized as follows: section 2 briefly describes the LIGO interferometers and the data from LIGO's S5 run. Section 3 defines the waveforms we seek and the associated assumptions we have made. In section 4 we briefly review the Hough-transform method. Section 5 describes the χ² test implemented for the analysis of the full S5 data. Section 6 gives a detailed description of the search pipeline and results. Upper limit computations are provided in section 7. The study of some features related to the χ²-veto is presented in section 8. Section 9 discusses variations, further improvements and capabilities of alternative searches. Section 10 concludes with a summary of the results.

During LIGO's S5 run the LIGO detector network consisted of a 4-km interferometer in Livingston, Louisiana (called L1) and two interferometers in Hanford, Washington, one a 4-km and another 2-km (H1 and H2, respectively). The S5 run spanned a nearly two-year period of data acquisition. This run started at 16:00 UTC on November 4, 2005 at Hanford and at 16:00 UTC on November 14, 2005 at Livingston Observatory; the run ended at 00:00 UTC on October 1, 2007. During this run, all three LIGO detectors had displacement spectral amplitudes very near their design goals of 1.1 × 10⁻¹⁹ m ⋅ Hz^−1/2 in their most sensitive frequency band near 150 Hz for the 4-km detectors and, in terms of GW strain, the H2 interferometer was roughly a factor of two less sensitive than the other two over most of the relevant band.

The data were acquired and digitized at a rate of 16384 Hz. Data acquisition was periodically interrupted by disturbances such as seismic transients (natural or anthropogenic), reducing the net running time of the interferometers. In addition, there were 1–2 week commissioning breaks to repair equipment and address newly identified noise sources. The resulting duty factors for the interferometers, defined as the fraction of the total run time when the interferometer was locked (i.e., all the interferometer control servos operating in their linear regime) and in its low configuration, were approximately 69% for H1, 77% for H2, and 57% for L1 during the first eight months. A nearby construction project degraded the L1 duty factor significantly during this early period of the S5 run. By the end of the S5 run, the cumulative duty factors had improved to 78% for H1, 79% for H2, and 66% for L1.

In the paper the data from each of the three LIGO detectors is used to search for continuous GW signals. In table 1 we provide the reference GPS initial and final times for the data collected for each detector, together with the number of hours of data used for the analysis, where each data segment used was required to contain at least 30 min of continuous interferometer operation.

Spinning neutron stars may generate continuous GWs through a variety of mechanisms. Independently of the specific mechanism, the emitted signal is a quasi-periodic wave whose frequency changes slowly during the observation time due to energy loss through GW emission, and possibly other mechanisms. The form of the received signal at the detector is

$$\begin{equation} h(t)=F_+\left( t,\psi \right) h_+ \left( t \right) + F_\times \left( t,\psi \right) h_\times \left( t \right) \end{equation} \tag{ 1 }$$

where t is time in the detector frame, ψ is the polarization angle of the wave and F_{+, ×} characterize the detector responses for the two orthogonal polarizations [29, 31]. For an isolated quadrupolar GW emitter, characterized by a rotating triaxial-ellipsoid mass distribution, the individual components h_{+, ×} have the form

$$\begin{equation} h_+ = h_0 \frac{1+\cos ^2\iota }{2} \cos \Phi (t) \quad \textrm {and} \quad h_\times = h_0 \cos \iota \sin \Phi (t), \end{equation} \tag{ 2 }$$

where ι describes the inclination of the source's rotation axis to the line of sight, h₀ is the wave amplitude and Φ(t) is the phase evolution of the signal.

For such a star, the GW frequency, f, is twice the rotation frequency and the amplitude h₀ is given by

$$\begin{equation} h_0 = \frac{4\pi ^2G}{c^4}\frac{I_{zz}f^2\varepsilon }{d}, \end{equation} \tag{ 3 }$$

where d is the distance to the star, I_zz is the principal moment of inertia with respect to its spin axis, ε the equatorial ellipticity of the star, G is Newton's constant and c is the speed of light.

Note that the search method used in this paper is sensitive to periodic signals from any type of isolated GW source, though we present upper limits in terms of h₀. Because we use the Hough method, only the instantaneous signal frequency in the detector frame, 2πf(t) = dΦ(t)/dt, needs to be calculated. This is given, to a very good approximation, by the non-relativistic Doppler expression:

$$\begin{equation} f(t)-\hat{f}(t) = \hat{f}(t) \frac{{\bf v}(t)\cdot {\bf n}}{c}, \end{equation} \tag{ 4 }$$

where $\hat{f}(t)$ is the instantaneous signal frequency in the solar system barycenter (SSB), v(t) is the detector velocity with respect to the SSB frame and n is the unit-vector corresponding to the sky location of the source. In this analysis, we search for $\hat{f}(t)$ signals well described by a nominal frequency f₀ at the start time of the S5 run t₀ and a constant first time derivative $\dot{f}$, such that

$$\begin{equation} \hat{f}(t) = f_0+ \dot{f} (t-t_0). \end{equation} \tag{ 5 }$$

These equations ignore corrections to the time interval t − t₀ at the detector compared with that at the SSB and relativistic corrections. These corrections are negligible for the search described here.

The Hough transform is a well known method for pattern recognition that has been applied to the search for continuous GWs. In this case the Hough transform is used to find hypothetical signals whose time-frequency evolution fits the pattern produced by the Doppler modulation of the detected frequency, due to the Earth's rotational and orbital motion with respect to the SSB, and the time derivative of the frequency intrinsic to the source. Further details can be found in [7, 21, 25, 26]; here we only give a brief summary.

The starting point for the Hough transform are N short Fourier transforms (SFTs). Each of these SFTs is digitized by setting a threshold ρ_th on the normalized power

$$\begin{equation} \rho _k = \frac{2|\tilde{x}_k|^2}{T_{{{{\rm coh}}}}S_n(f_k)} \,. \end{equation} \tag{ 6 }$$

Here ${\tilde{x}_k}$ is the discrete Fourier transform of the data, the frequency index k corresponds to a physical frequency of f_k = k/T_coh, S_n(f_k) is the single sided power spectral density of the detector noise and T_coh is the time baseline of the SFT. The kth frequency bin is selected if ρ_k ⩾ ρ_th, and rejected otherwise. In this way, each SFT is replaced by a collection of zeros and ones called a peak-gram. This is the simplest method of selecting frequency bins, for which the optimal choice of the threshold ρ_th is 1.6 [21]. Alternative conditions could be imposed [32–34], that might be more robust against spectral disturbances.

For our choice, the probability that a frequency bin is selected is $q = {\rm e}^{-\rho _{{{\rm th}}}}$ for Gaussian noise and η, given by

$$\begin{equation} \eta = q\left\lbrace 1+\frac{\rho _{{{\rm th}}}}{2}\lambda _k + \mathcal {O}(\lambda _k^2) \right\rbrace \end{equation} \tag{ 7 }$$

is the corresponding probability in the presence of a signal. λ_k is the signal to noise ratio (SNR) within a single SFT, and for the case when there is no mismatch between the signal and the template:

$$\begin{equation} \lambda _k = \frac{4|\tilde{h}(f_k)|^2}{T_{{{{\rm coh}}}}S_n(f_k)} \end{equation} \tag{ 8 }$$

with $\tilde{h}(f)$ being the Fourier transform of the signal h(t).

Several flavors of the Hough transform have been developed [21, 25, 35] and used for different searches [7, 10, 14]. The Hough transform is used to map points from the time-frequency plane of our data (understood as a sequence of peak-grams) into the space of the source parameters. Each point in parameter space corresponds to a pattern in the time-frequency plane, and the Hough number count n is the weighted sum of the ones and zeros, $n_k^{(i)}$, of the different peak-grams along this curve. For the 'weighted Hough' this sum is computed as

$$\begin{equation} n = \sum _{i=0}^{N-1} w^{(i)}_{k} n^{(i)}_{k}, \end{equation} \tag{ 9 }$$

where the choice of weights is optimal, in the sense of [25], if defined as

$$\begin{equation} w^{(i)}_{k} \propto \frac{1}{S^{(i)}_{k}}\big\lbrace \big(F_{+1/2}^{(i)}\big)^2 + \big(F_{\times 1/2}^{(i)}\big)^2\big\rbrace , \end{equation} \tag{ 10 }$$

where $F_{+1/2}^{(i)}$ and $F_{\times 1/2}^{(i)}$ are the values of the beam pattern functions at the mid point of the ith SFT and are normalized according to

$$\begin{equation} \sum _{i=0}^{N-1} w^{(i)}_{k} = N\,. \end{equation} \tag{ 11 }$$

The natural detection statistic is the significance (or critical ratio) defined as:

$$\begin{equation} s= \frac{n-\langle n\rangle }{\sigma } , \end{equation} \tag{ 12 }$$

where 〈n〉 and σ are the expected mean and standard deviation for pure noise. Furthermore, the relation between the significance and the false alarm probability α, in the Gaussian approximation [21], is given by

$$\begin{equation} s_{{{\rm th}}}= \sqrt{2}\textrm {erfc}^{-1}(2\alpha )\,. \end{equation} \tag{ 13 }$$

χ² time-frequency discriminators are commonly used for GW detection. Originally, they were designed for broadband signals with a known waveform in a data stream [36]. But they can be adapted for narrowband continuous signals, as those expected from rapidly rotating neutron stars. The essence of these tests is to 'break up' the data (in time or frequency domain) and to see if the response in each chunk is consistent with what would be expected from the purported signal.

In this paper, a chi-square test is implemented as a veto, in order to reduce the number of candidates in the analysis of the full S5 data. The idea for this χ² discriminator is to split the data into p non-overlapping chunks, each of them containing a certain number of SFTs {N₁, N₂, ..., N_p}, such that

$$\begin{equation} \sum _{j=1}^p N_j =N , \end{equation} \tag{ 14 }$$

and analyze them separately, obtaining the Hough number-count n_j which, for the same pattern across the different chunks, would then satisfy

$$\begin{equation} \sum _{j=1}^p n_j =n , \end{equation} \tag{ 15 }$$

where n is the total number-count for a given point in parameter space. The χ² statistic will look along the different chunks to see if the number count accumulates in a way that is consistent with the properties of the signal and the detector noise. Small values of χ² are consistent with the hypothesis that the observed significance arose from a detector output which was a linear combination of Gaussian noise and the continuous wave signal. Large values of χ² indicate either the signal did not match the template or that the detector noise was non-Gaussian.

In the following subsections we derive a χ² discriminator for the different implementations of the Hough transform and show how the veto curve was derived for LIGO S5 data.

5.1. The standard Hough

In the simplest case in which all weights are set to unity, the expected value and variance of the number count are

$$\begin{equation} \langle n\rangle =N\eta , \quad \sigma ^2_n =N\eta (1-\eta ) , \end{equation} \tag{ 16 }$$

$$\begin{equation} \langle n_j\rangle = N_j\eta = N_j\frac{\langle n\rangle }{N} , \quad \sigma ^2_{n_j} =N_j\eta (1-\eta ) . \end{equation} \tag{ 17 }$$

Consider the p quantities defined by

$$\begin{equation} \Delta n_j\equiv n_j-\frac{N_j}{N}n . \end{equation} \tag{ 18 }$$

With this definition, it holds true that

$$\begin{equation} \langle \Delta n_j\rangle =0 , \quad \sum _{j=1}^p \Delta n_j=0 , \quad \langle n_j n\rangle = \frac{N_j}{N} \langle n^2\rangle , \end{equation} \tag{ 19 }$$

and the expectation value of the square of Δn_j is

$$\begin{equation} \langle (\Delta n_j)^2\rangle =\left( 1-\frac{N_j}{N}\right) N_j\eta (1-\eta ) . \end{equation} \tag{ 20 }$$

Therefore we can define the χ² discriminator statistic by

$$\begin{equation} \chi ^2(n_1, \ldots , n_p) = \sum _{j=1}^p {\frac{(\Delta n_j)^2}{\sigma ^2_{n_j}}}= \sum _{j=1}^p {\frac{(n_j-nN_j/N)^2}{N_j\eta (1-\eta )} } . \end{equation} \tag{ 21 }$$

This corresponds to a χ²-distribution with p − 1 degrees of freedom. To implement this discriminator, we need to measure, for each point in parameter space, the total number-count n, the partial number-counts n_j and assume a constant value of η = n/N.

5.2. The weighted Hough

In the case of the weighted Hough the result given by equation (21) can be generalized. Let I_j be the set of SFT indices for each different p chunks, thus the mean and variance of the number-count become

$$\begin{equation} \langle n_j \rangle = \sum _{i\in I_j} w_i\eta _i \quad \langle n \rangle = \sum _{j=1}^p \langle n_j \rangle \quad \sigma ^2_{n_j} = \sum _{i\in I_j} w_i^2\eta _i(1-\eta _i) \end{equation} \tag{ 22 }$$

and we can define

$$\begin{equation} \Delta n_j\equiv n_j-n\frac{\sum _{i\in I_j} w_i\eta _i}{\sum _{i=1}^N w_i\eta _i} , \end{equation} \tag{ 23 }$$

so that 〈Δn_j〉 = 0, $\sum _{j=1}^p \Delta n_j=0$. Hence, the χ² discriminator would now be:

$$\begin{equation} \chi ^2 = \sum _{j=1}^p {\frac{(\Delta n_j)^2}{\sigma ^2_{n_j}} } = \sum _{j=1}^p {\frac{\big(n_j-n \big(\sum _{i\in I_j} w_i\eta _i\big)\big/\big(\sum _{i=1}^N w_i\eta _i\big)\big)^2}{\sum _{i\in I_j} w_i^2\eta _i(1-\eta _i)} } . \end{equation} \tag{ 24 }$$

In a given search, we can compute the $\sum _{i\in I_j} w_i$, $\sum _{i\in I_j} w_i^2$ for each of the p chunks, but the different η_i values cannot be measured from the data itself because they depend on the exact SNR for each single SFT as defined in equations (7) and (8). For this reason, the discriminator we proposed is constructed by replacing η_i → η*, where η* = n/N. In this way, from equation (24) we get

$$\begin{equation} \chi ^2 \approx \sum _{j=1}^p {\frac{\big(n_j-n \big(\sum _{i\in I_j} w_i\big)/N\big)^2}{\eta ^*(1-\eta ^*) \sum _{i\in I_j} w_i^2 } } . \end{equation} \tag{ 25 }$$

In principle, one is free to choose the different p chunks of data as one prefers, but it is reasonable to split the data into segments in such a way that they would contribute a similar relative contribution to the total number count. Therefore we split the SFT data in such a way that the sum of the weights in each block satisfies

$$\begin{equation} \sum _{i\in I_j} \omega _{i} \approx \frac{N}{p} . \end{equation} \tag{ 26 }$$

Further details and applications of this χ² on LIGO S4 data can be found in [27].

5.3. The S5 χ² veto curve

We study the behavior of this χ² discriminator in order to characterize the χ²-significance plane in the presence of signals and derive empirically the veto curve.

For this purpose we use the full LIGO S5 SFT data, split into both years, in the same way that is done in the analysis, and inject a large number of Monte Carlo simulated continuous GW signals into the data, varying the amplitude, frequency, frequency derivative, sky location, as well as the nuisance parameters cos ι, ψ and ϕ₀ of the signals. Those injections are analyzed with the multi-interferometer Hough code using the same grid resolution in parameter space as is used in the search.

To characterize the veto curve, nine 0.25 Hz bands, spread in frequency and free of known large spectral disturbances have been selected. These are: 102.5, 151, 190, 252.25, 314.1, 448.5, 504.1, 610.25 and 710.25 Hz. Monte Carlo injections in those bands have been performed separately in the data from both years. Since the results were comparable for both years a single veto curved is derived.

In total 177 834 injections are considered with a significance value lower than 70. The results of these injections in terms of (s, χ²) are presented in figure 1. The χ² values obtained correspond to those by splitting the data in p = 16 segments.

Zoom In Zoom Out Reset image size

Then we proceed as follows: first we sort the points with respect to the significance, and we group them in sets containing 1000 points. For each set we compute the mean value of the significance, the mean of the χ² and its standard deviation. With these reduced set of points we fit two power laws p − 1 + a s^c and $\sqrt{2p-2}+b \, s^d$ to the (mean s, mean χ²) and (mean s, std χ²) respectively, obtaining the following coefficients (with 95% confidence bounds):

$$\begin{eqnarray*} a = 0.3123 & \quad & (0.305, 0.3195) \nonumber \\ c = 1.777 & \quad & (1.77, 1.783) \nonumber \\ b = 0.1713 & \quad & (0.1637, 0.1789) \nonumber \\ d = 1.621 & \quad & (1.609, 1.633). \nonumber \end{eqnarray*}$$

The veto curve we will use in this analysis corresponds to the mean curve plus five times the standard deviation

$$\begin{equation} \bar{\chi }^2= p-1 + 0.3123 s^{1.777} + 5 (\sqrt{2p-2} + 0.1713 s^{1.621}). \end{equation} \tag{ 27 }$$

This curve vetoes 25 of the 177 834 injections considered with significance lower than 70, that could translate into a false dismissal rate of 0.014. In figure 1 we show the fitted curves and the $\bar{\chi }^2$ veto curve compared to the result of the injections.

In this paper, we use a new pipeline to analyze the data from the S5 run of the LIGO detectors to search for evidence of continuous GWs, that might be radiated by nearby unknown rapidly spinning isolated neutron stars. Data from each of the three LIGO interferometers is used to perform the all-sky search. The key difference from previous searches is that, starting from 30 min SFTs, we perform a multi-interferometer search analyzing separately the two years of the S5 run, and we study coincidences among the source candidates produced by the first and second years of data. Furthermore, we use a χ² test adapted to the Hough transform searches to veto potential candidates. The pipeline is shown schematically in figure 2.

Zoom In Zoom Out Reset image size

A separate search was run for each successive 0.25 Hz band within the frequency range 50–1000 Hz and covering frequency time derivatives in the range −8.9 × 10⁻¹⁰ Hz s⁻¹ to zero. We use a uniform grid spacing equal to the size of a SFT frequency bin, δf = 1/T_coh = 5.556 × 10⁻⁴ Hz. The resolution $\delta \dot{f}$ is given by the smallest value of $\dot{f}$ for which the intrinsic signal frequency does not drift by more than a frequency bin during the observation time T_obs in the first year: $\delta \dot{f} = \delta f / T_{{{{\rm obs}}}}\sim 1.8\times 10^{-11} \,\mathrm{Hz} \mathrm{\hbox{ }s}^{-1}$. This yields 51 spin-down values for each frequency. $\delta \dot{f}$ is fixed to the same value for the search on the first and the second year of S5 data. The sky resolution, δθ, is frequency dependent, as given by equation (4.14) of [21], that we increase by a factor 2. As explained in detail in section V.B.1 of [10], the sky-grid spacing can be increased with a negligible loss in SNR, and for previous PowerFlux searches [10, 13, 17] a factor 5 of increase was used in some frequency ranges to analyze LIGO S4 and S5 data.

The set of SFTs are generated directly from the calibrated data stream, using 30-min intervals of data for which the interferometer is operating in what is known as science mode. With this requirement, we search 32295 SFTs from the first year of S5 (11402 from H1, 12195 from H2 and 8698 from L1) and 35401 SFTs from the second year (12590 from H1, 12178 from H2 and 10633 from L1).

6.1. A two-step hierarchical Hough search

The approach used to analyze each year of data is based on a two-step hierarchical search for continuous signals from isolated neutron stars. In both steps, the weighted Hough transform is used to find signals whose frequency evolution fits the pattern produced by the Doppler shift and the spin-down in the time-frequency plane of the data. The search is done by splitting the frequency range in 0.25 Hz bands and using the SFTs from multiple interferometers.

In the first stage, and for each 0.25 Hz band, we break up the sky into smaller patches with frequency dependent size in order to use the look up table approach to compute the Hough transform, which greatly reduces the computational cost. The look up table approach benefits from the fact that, according to the Doppler expression (4), the set of sky positions consistent with a given frequency bin f_k at a given time correspond to annuli on the celestial sphere centered on the velocity vector v(t). In the look up table approach, we precompute all the annuli for a given time and a given search frequency mapped on the sky search grid. Moreover, it turns out that the mapped annuli are relatively insensitive to changes in frequency and can therefore be reused a large number of times. The Hough map is then constructed by selecting the appropriate annuli out of all the ones that have been found and adding them using the corresponding weights. A detailed description of the look up table approach with further details of implementation choices can be found in [21].

But limitations on the memory of the computers constrain the volume of data (i.e., the number of SFTs) that can be analyzed at once and the parameter space (e.g., size and resolution of the sky-patches and number of spin-down values) we can search over. For this reason, in this first stage, we select the best 15000 SFTs (according to the noise floor and the beam pattern functions) for each frequency band and each sky-patch and apply the Hough transform on the selected data. The size of the sky-patches ranges from ∼0.4 rad × 0.4 rad at 50 Hz to ∼0.07 rad × 0.07 rad at 1 kHz and we calculate the weights only for the center of each sky-patch. This was set in order to ensure that the memory usage will never exceed the 0.8 GB and this search could run on the Merlin/Morgane dual compute cluster at the Albert Einstein Institute¹²⁸. A top-list keeping the best 1000 candidates is produced for each 0.25 Hz band for the all-sky search.

Figure 3 shows the histograms of the percentage of SFTs that each detector contributes for the different sky locations for a band at 420 Hz for the first year of S5 data. At this particular frequency, the detector that contributes the most is H1 between 44–64%, giving the maximum contribution near the poles, L1 contributes between 28.1–45.5% with its maximum around the equator, and H2 contributes at most 21.7% of the SFTs. As shown in figure 1 in [37], the maximum contribution of H2 corresponds to those sky regions where L1 contributes the least. If SFT selection had been based only upon the weights due to the noise floor, the H2 detector would not have contributed at all in this first stage.

Zoom In Zoom Out Reset image size

In a second stage, we compute the χ² value for all the candidates in the top-list in each 0.25 Hz band. This is done by dividing the data into 16 chunks and summing weighted binary zeros or ones along the expected path of the frequency evolution of a hypothetical periodic GW signal in the digitized time-frequency plane of our data. Since there are no computational limitations, we use the complete set of available SFTs from all three interferometers, and we also get a new value of the significance using all the data. In this way we reduce the mismatch of the template, since the number count is obtained without the roundings introduced by the look up table approach and the weights are computed for the precise sky location and not for the center of the corresponding patch. All these refinements contribute also to a potential improvement of sensitivity when a threshold is subsequently applied to the recomputed significance (described below).

Figure 4 shows the maximum-significance value in each 0.25 Hz band obtained for the first and second years of S5 data.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

6.2. The post-processing

After the multi-interferometer Hough search is performed on each year of S5 data between 50 and 1000 Hz, a top list keeping the best 1000 candidates is produced for each 0.25 Hz band. This step yields 3.8 × 10⁶ candidates for each year. The post-processing of these results has the following steps:

• (i)  
  Remove those 0.25 Hz bands that are affected by power lines or violin modes.A total of 96 bands are removed. These bands are given in table 2.
• (ii)  
  Remove all the 0.25 Hz bands for which the χ² vetoes more than a 95% of the elements in the top list.Figure 5 shows the χ² veto level for all the frequency bands for both years. With this criterion, 144 and 131 0.25 Hz bands would be vetoed for the first and second year of data respectively. These first two steps leave a total 3548 bands in which we search for coincidence candidates and set upper limits; a total of 252 bands were discarded.
• (iii)  
  Set a threshold on the significance.Given the relation of the Hough significance and the Hough false alarm probability (see equation (13)), we set a threshold on the candidate's significance that corresponds to a false alarm of 1/(number of templates) for each 0.25 Hz band. Figure 6 shows the value of this threshold at different frequencies.
• (iv)  
  Apply the χ² veto.From the initial 3.8 × 10⁶ elements in the top list for each year, after excluding the noisy bands, applying the χ² veto and setting a threshold on the significance, the number of candidates remaining in the 3548 'clean' bands are 31 427 for the first year and 50 832 for the second year. Those are shown in figure 7.
• (v)  
  Selection of coincident candidates.For each of the four parameters: frequency, spin-down and sky location, we set the coincidence window with a size equal to five times the grid spacing used in the search and centered on the values of the candidates parameters. Therefore the coincidence window always contains 625 cells in parameter space, with frequency-dependent size, according to the search grid. This window is computed for each of the candidates selected from the first year of data and then we look for coincidences among the candidates of the second year, making sure to translate their frequency to the reference time of the starting time of the run, taking into account their spin-down values. Extensive analysis of software injected signals, in different frequency bands, have been used to determine the size of this coincidence window. This was done by comparing the parameters of the most significant candidates of the search, using the same pipeline, in both years of data.With this procedure, we obtain 135 728 coincidence pairs, corresponding to 5823 different candidates of the first year that have coincidences with 7234 different ones of the second year. Those are displayed in figure 8. All those candidates cluster in frequency in 34 groups. The most significant outlier at 108.857 Hz corresponds to a simulated pulsar signal injected into the instrument as a test signal. The most significant events in each cluster are shown in tables 3 and 4. Notice how with this coincidence step the overall number of candidates has been reduced by a factor 5.4 for the first year and a factor 7.0 for the second year. Furthermore, without the coincidence step, the candidates are spread over all frequencies, whereas the surviving coincident candidates are clustered in a few small regions, illustrating the power of this procedure on real data.

Table 2. Initial frequency of the 0.25 Hz bands excluded from the search.

Excluded Bands (Hz)	Description
[n60 − 0.25, n60 + 0.25] n = 1 to 16	Power lines
[343.0, 344.75]	Violin modes
[346.5, 347.75]	Violin modes
[348.75, 349.25]	Violin modes
[685.75, 689.75]	Violin mode harmonics
[693.0, 695.5]	Violin mode harmonics
[697.5, 698.75]	Violin mode harmonics

Table 3. Summary of first year coincidence candidates, including the frequency band, the number of candidates in each cluster and showing the details of the most significant candidate in each of the 34 clusters. Shown are the significance s, the χ² value, the detected frequency at the start of the run (SSB frame) f₀, the spin-down $\dot{f}$, and the sky position (RA, dec).

	Band (Hz)	Num.	s	χ2	f0 (Hz)	$\dot{f}$ (Hz s−1)	RA (rad)	dec (rad)
1	50.001–50.003	12	7.103	64.115	50.0022	0	−1.84	0.69
2	50.997–51.004	7	5.431	57.322	51.0028	0	−2.73	0.94
3	52.000–52.016	44	6.886	54.108	52.0139	−21.4e-11	−1.99	−0.08
4	52.786–52.793	6	6.094	40.089	52.7911	0	0	1.33
5	53.996–54.011	1136	11.880	60.583	54.0039	−10.7e-11	2.68	−1.21
6	54.996–55.011	82	6.995	48.289	55.0067	0	3.00	−0.47
7	55.749–55.749	7	5.940	37.902	55.7489	0	−0.23	1.12
8	56.000–56.016	167	8.085	66.130	56.0056	−7.1e-11	0.20	0.22
9	56.997–57.011	1370	24.751	175.459	57.0028	−12.5e-11	−2.10	1.40
10	58.000–58.015	89	7.901	61.360	58.0128	−35.7e-11	−2.32	0.43
11	61.994–62.000	195	6.597	26.445	61.9989	0	0.98	0.21
12	62.996–63.009	1324	23.188	131.619	63.0028	−8.9e-11	0	−1.40
13	64.996–65.000	101	6.252	52.374	64.9994	−5.3e-11	−0.17	0.03
14	65.378–65.381	3	5.586	29.000	65.3806	−16.1e-11	−2.17	1.15
15	65.994–66.013	919	11.662	35.587	65.9994	−7.1e-11	0.15	1.38
16	67.006–67.006	1	5.801	11.950	67.0056	−17.8e-11	−1.13	−0.20
17	67.993–68.009	39	6.061	44.319	68.0017	−14.3e-11	−1.55	0.66
18	72.000–72.000	4	5.604	17.668	72.0000	−1.8e-11	1.57	−1.11
19	86.002–86.024	14	6.786	47.054	86.0150	−17.8e-11	1.78	1.09
20	90.000–90.000	2	5.554	58.338	90.0000	−3.6e-11	1.54	−1.05
21	108.857–108.860	50	64.850	552.161	108.8570	0	3.10	−0.60
22	111.998–111.998	1	5.673	46.493	111.9980	−7.1e-11	−0.54	1.18
23	118.589–118.613	18	7.072	52.181	118.5990	−57.1e-11	2.90	−0.46
24	160.000–160.000	1	5.571	18.847	160.0000	0	1.56	−1.15
25	178.983–179.026	21	7.483	43.380	179.0010	−3.6e-11	−1.59	1.17
26	181.000–181.038	8	6.309	13.174	181.0170	−8.9e-11	−1.12	−0.89
27	192.000–192.002	5	7.976	41.042	192.0000	−1.8e-11	−1.51	1.17
28	341.763–341.765	3	6.292	39.340	341.7630	−1.8e-11	−1.63	1.21
29	342.680–342.684	5	6.113	36.893	342.6800	−3.6e-11	1.47	−1.14
30	345.721–345.724	17	6.835	38.540	345.7230	−12.5e-11	−1.07	1.44
31	346.306–346.316	9	6.973	30.298	346.3070	−3.6e-11	1.32	−1.13
32	394.099–394.100	3	10.617	95.063	394.1000	0	−1.58	1.16
33	575.163–575.167	57	26.058	146.929	575.1640	−1.8e-11	−2.53	0.06
34	671.728–671.733	101	13.878	132.197	671.7290	0	1.54	−1.17

Table 4. Summary of 2nd year coincidence candidates, showing the details of the most significant candidate in each of the 34 clusters.

	Band (Hz)	Num.	s	χ2	f0 (Hz)	$\dot{f}$ (Hz s−1)	RA (rad)	dec (rad)
1	50.001–50.004	10	9.345	77.005	50.0006	−5.3e-11	1.25	−1.05
2	50.993–51.003	16	7.106	17.670	51.0011	−1.8e-11	−1.95	1.32
3	52.000–52.012	38	8.800	30.888	52.0094	−17.8e-11	−2.17	−0.08
4	52.784–52.792	14	13.269	128.801	52.7911	−7.1e-11	0.77	1.31
5	53.996–54.006	1380	23.431	69.473	54.0033	−14.3e-11	−1.92	1.31
6	54.995–55.008	405	11.554	73.521	55.0056	−5.3e-11	−2.98	0.13
7	55.748–55.749	46	7.102	34.151	55.7489	−1.8e-11	−0.10	0.13
8	56.000–56.008	161	11.478	64.077	56.0006	−1.8e-11	1.58	−1.13
9	56.996–57.004	1478	47.422	412.817	56.9989	0	−0.31	1.43
10	58.000–58.008	166	11.675	65.764	58.0000	0	1.59	−1.13
11	61.993–62.001	400	9.361	56.108	61.9989	0	−1.16	1.02
12	62.997–63.004	1138	42.432	360.615	63.0039	−5.3e-11	−1.51	1.41
13	64.994–65.000	288	12.557	119.361	64.9978	0	0.61	−1.14
14	65.377–65.378	2	5.777	34.580	65.3772	−14.3e-11	−1.84	1.03
15	65.995–66.006	1162	20.864	67.221	66.0022	−5.3e-11	−1.51	1.40
16	66.999–66.999	5	6.207	19.546	66.9989	−17.8e-11	−1.08	0.03
17	67.994–68.007	199	11.682	82.662	68.0006	−1.8e-11	1.65	−1.14
18	71.999–72.000	10	6.802	45.727	72.0000	−1.8e-11	1.41	−1.14
19	86.002–86.013	15	7.368	54.707	86.0094	−10.7e-11	2.10	0.75
20	89.999–90.000	4	9.008	86.486	90.0000	0	1.54	−1.19
21	108.857–108.858	27	77.157	1613.230	108.8580	−5.3e-11	2.99	−0.71
22	111.996–111.996	1	5.775	59.202	111.9960	−7.1e-11	−0.37	1.00
23	118.579–118.589	19	7.398	69.372	118.5820	−41.0e-11	2.79	−0.68
24	160.000–160.000	2	7.898	15.079	160.0000	0	1.60	−1.17
25	178.984–179.014	24	7.089	33.229	179.0010	−10.7e-11	−1.79	1.35
26	180.998–181.019	8	6.587	28.605	181.0180	−28.5e-11	−0.39	1.07
27	191.999–192.001	8	6.582	45.301	192.0000	−3.6e-11	1.51	−1.19
28	341.762–341.764	19	7.859	71.736	341.7630	−3.6e-11	−1.63	1.17
29	342.677–342.680	19	7.569	45.255	342.6790	−14.3e-11	1.47	−1.11
30	345.718–345.721	12	7.890	44.570	345.7200	−12.5e-11	−0.91	1.38
31	346.303–346.309	14	7.803	45.756	346.3070	0	1.46	−1.21
32	394.099–394.100	2	9.877	95.232	394.1000	0	−1.58	1.16
33	575.163–575.165	53	40.576	415.830	575.1640	−1.8e-11	−2.53	0.06
34	671.729–671.732	87	13.600	116.261	671.7320	−1.8e-11	1.59	−1.15

Noise lines were identified by previously performed searches [13, 14, 17, 18, 38] as well as the search described in this paper. Several techniques were used to identify the causes of outliers, including the calculation of the coherence between the interferometer output channel and physical environment monitoring channels and the computation of high resolution spectra. A dedicated analysis code 'FScan' [39] was also created specifically for identification of instrumental artifacts. Problematic noise lines were recorded and monitored throughout S5.

In addition, a number of particular checks were performed on the coincidence outliers, including: a detailed study of the full top-list results, for those 0.25 Hz bands where the candidates were found—in order to check if candidates are more dominant in a given year, or if they cluster in certain regions of parameter space; and a second search using the data of the two most sensitive detectors, H1 and L1 separately—in order to see if artifacts could be associated to a given detector, consistent with the observed spectra.

All of the 34 outliers were investigated and were all traced to instrumental artifacts or hardware injections (see details in table 5). Hence the search did not reveal any true continuous GW signals.

Table 5. Description of the coincidence outliers, together with the maximum value of the significance in both years.

	Bands (Hz)	s 1y	s 2y	Comment
1	50.001–50.004	7.103	9.345	L1 1 Hz Harmonic from control/data acquisition system
2	50.993–51.004	5.431	7.106	L1 1 Hz Harmonic from control/data acquisition system
3	52.000–52.016	6.886	8.800	L1 1 Hz Harmonic from control/data acquisition system
4	52.784–52.793	6.094	13.269	Instrumental line in H1
5	53.996–54.011	11.880	23.431	Pulsed heating sideband on 60 Hz mains
6	54.995–55.011	6.995	11.554	1 Hz Harmonic from control/data acquisition system
7	55.748–55.749	5.940	7.102	Instrumental line in L1
8	56.000–56.016	8.085	11.478	L1 1 Hz Harmonic from control/data acquisition system
9	56.996–57.011	24.751	47.422	Pulsed heating sideband on 60 Hz mains
10	58.000–58.015	7.901	11.675	L1 1 Hz Harmonic from control/data acquisition system
11	61.993–62.001	6.597	9.361	L1 1 Hz Harmonic from control/data acquisition system
12	62.996–63.009	23.188	42.432	Pulsed heating sideband on 60 Hz mains
13	64.994–65.000	6.252	12.557	L1 1 Hz Harmonic from control/data acquisition system
14	65.377–65.381	5.586	5.777	Instrumental line in L1—member of offset 1 Hz comb
15	65.994–66.013	11.662	20.864	Pulsed sideband on 60 Hz mains
16	66.999–67.006	5.801	6.207	L1 1 Hz Harmonic from control/data acquisition system
17	67.993–68.009	6.061	11.682	1 Hz Harmonic from control/data acquisition system
18	71.999–72.000	5.604	6.802	1 Hz Harmonic from control/data acquisition system
19	86.002–86.024	6.786	7.368	Instrumental line in H1
20	89.999–90.000	5.554	9.008	Instrumental line in H1
21	108.857–108.860	64.850	77.157	Hardware injection of simulated signal (ip3)
22	111.996–111.998	5.673	5.775	16 Hz harmonic from data acquisition system
23	118.579–118.613	7.072	7.398	Sideband of mains at 120 Hz
24	160.000–160.000	5.571	7.898	16 Hz harmonic from data acquisition system
25	178.983–179.026	7.483	7.089	Sideband of mains at 180 Hz
26	180.998–181.038	6.309	6.587	Sideband of mains at 180 Hz
27	191.999–192.002	7.976	6.582	16 Hz harmonic from data acquisition system
28	341.762–341.765	6.292	7.859	Sideband of suspension wire resonance in H1
29	342.677–342.684	6.113	7.569	Sideband of suspension wire resonance in H1
30	345.718–345.724	6.835	7.890	Sideband of suspension wire resonance in H1
31	346.303–346.316	6.973	7.803	Sideband of suspension wire resonance in H1
32	394.099–394.100	10.617	9.877	Sideband of calibration line at 393.1 Hz in H1
33	575.163–575.167	26.058	40.576	Hardware injection of simulated signal (ip2)
34	671.728–671.733	13.878	13.600	Instrumental line in H1

The analysis of the Hough search presented here has not identified any convincing continuous GW signal. Hence, we proceed to set upper limits on the maximum intrinsic GW strain h₀ that is consistent with our observations for a population of signals described by an isolated triaxial rotating neutron star.

As in the previous S2 and S4 searches [7, 10], we set a population-based frequentist upper limit, assuming random positions in the sky, in the GW frequency range [50, 1000] Hz and with spin-down values in the range −8.9 × 10⁻¹⁰ Hz s⁻¹ to zero. The rest of the nuisance parameters, cos ι, ψ and ϕ₀, are assumed to be uniformly distributed. As commonly done in all-sky, all-frequency searches, the upper limits are given in different frequency sub-bands, here chosen to be 0.25 Hz wide. Each upper limit is based on the most significant event from each year in its 0.25 Hz band. Our goal is to find the value of h₀ (denoted $h_0^{95\%}$) such that 95% of the signal injections at this amplitude would be recovered by our search and are more significant than the most significant candidate from the actual search in that band, thus providing the 95% confidence all-sky upper limit on h₀.

Our procedure for setting upper-limits uses partial Monte Carlo signal injection studies, using the same search pipeline as described above, together with an analytical sensitivity estimation. As in the previous S4 Hough search [10], upper limits can be computed accurately without extensive Monte Carlo simulations. Up to a constant factor C, that depends on the grid resolution in parameter space, they are given by

$$\begin{equation} h_0^{95\%} = C \left( \frac{1}{ \sum _{i=0}^{N-1}(S^{}_{i})^{-2}}\right)^{1/4}\sqrt{\frac{ \mathcal{S}}{T_{{{{\rm coh}}}}}} \,. \end{equation} \tag{ 28 }$$

where

$$\begin{equation} \mathcal{S}= \textrm {erfc}^{-1}(2\alpha _{{{\rm H}}}) + \textrm {erfc}^{-1}(2\beta _{{{\rm H}}}) , \end{equation} \tag{ 29 }$$

S_i is the average value of the single sided power spectral noise density of the ith SFT in the corresponding frequency sub-band, α_H is the false alarm and β_H the false dismissal probability.

The utility of this fit is that having determined the value of C in a small frequency range, it can be extrapolated to cover the full bandwidth without performing any further Monte Carlo simulations. Figure 9 shows the value of the constant C for a number of 0.1 Hz frequency bands. More precisely, this is the ratio of the upper limits measured by means of Monte-Carlo injections in the multi-interferometer Hough search to the quantity $h_0^{95\%}/C$ as defined in equation (28). The value of $\mathcal{S}$ is computed using the false alarm α_H corresponding to the observed loudest event, in a given frequency band, and a false dismissal rate β_H = 0.05, in correspondence to the desired confidence level of 95%, i.e., $\mathcal{S}\rightarrow s^*/ \sqrt{2} + \textrm {erfc}^{-1}(0.1)$, where s* is the highest significance value in the frequency band. This yields a scale factor C of 8.32 ± 0.19 for the first year and 8.25 ± 0.16 for the second year of S5. With these values we proceed to set the upper limits for all the frequency bands. The validity of equation (28) was studied in [10] using LIGO S4 data. In that paper upper limits were measured for each 0.25 Hz frequency band from 100 to 1000 Hz using Monte Carlo injections and compared with those prescribed by this analytical approximation. Such comparison study showed that the values obtained using equation (28) have an error smaller than 5% for bands free of large instrumental disturbances. For an in-depth study of how to analytically estimate the sensitivity of wide parameter searches for GW pulsars, we refer the reader to [30].

Zoom In Zoom Out Reset image size

The 95% confidence all-sky upper limits on h₀ from this multi-interferometer search for each year of S5 data are shown in figure 10. The best upper limits correspond to 1.0 × 10⁻²⁴ for the first year of S5 in the 158–158.25 Hz band, and 8.9 × 10⁻²⁵ for the second year in the 146.5–146.75 Hz band. There is an overall 15% calibration uncertainty on these upper limits. No upper limits are provided in the 252 vetoed bands, that were excluded from the coincidence analysis, since the analytical approximation would not be accurate enough. These excluded frequency bands are marked in the figure.

Zoom In Zoom Out Reset image size

Figure 11 provides the maximum astrophysical reach of our search for each year of the S5 run. The top panel shows the maximum distance to which we could have detected a source emitting a continuous wave signal with strain amplitude $h_0^{95\%}$. The bottom panel does not depend on any result from the search. It shows the corresponding ellipticity values as a function of frequency. For both plots the source is assumed to be spinning down at the maximum rate considered in the search −8.9 × 10⁻¹⁰ Hz s⁻¹, and emitting in GWs all the energy lost. This follows formulas in paper [10] and assumes the canonical value of 10³⁸ kg m² for I_zz in equation (3).

Zoom In Zoom Out Reset image size

Around the frequencies of greatest sensitivity, we are sensitive to objects as far away as 1.9 and 2.2 kpc for the first and second year of S5 and with an ellipticity ε around 10⁻⁴. Normal neutron stars are expected to have ε less than 10⁻⁵ [40, 41]. Such plausible value of ε could be detectable by a search like this if the object were emitting at 350 Hz and at a distance no further than 750 pc. For a source of fixed ellipticity and frequency, this search had a bit less range than the Einstein@Home search on the same data [18].

A novel feature of the search presented here is the implementation of the χ² veto. It is worth mentioning that this discriminator has been able to veto all the violin modes present in the data and many other narrow instrumental artifacts. Figure 12 demonstrates how well the χ² veto used works on those frequency bands affected by violin modes.

Zoom In Zoom Out Reset image size

As part of the testing and validation of search pipelines and analysis code, simulated signals are added into the interferometer length control system to produce mirror motions similar to what would be generated if a GW signal were present. These are the so-called hardware-injected pulsars. During the S5 run ten artificial pulsars were injected. Four of these pulsars: P2, P3, P5 and P8, at frequencies 575.16, 108.85, 52.81 and 193.4 Hz respectively, were strong enough to be detected by the multi-interferometer Hough search (see table III in [18] for the detailed parameters). The hardware injections were not active all the time, having a duty factor of about 60%.

The fact that these signals were not continuously present in the data caused the χ² test to veto most of the templates associated with them, since they did not behave like the signals we were looking for. In particular, for the second year of S5, the elements of the top-list in frequency band containing P8 were vetoed by the χ² test at the 99.4% level, and therefore that band was excluded from the analysis. The bands containing injected pulsars P2 and P3 were vetoed at the 87.7% and 94.5% level respectively, including the most significant events. Figure 13 shows the behavior of the χ² veto for the 0.25 Hz band starting at 108.75 Hz that contains pulsar P3.

Zoom In Zoom Out Reset image size

In the frequency band 52.75–53.0 Hz, the candidates in the top-list were all produced by the 52.79 Hz instrumental artifact present in H1 and consequently the search failed to detect P5. This suggests that, in future analysis, smaller frequency intervals should be used to produce the top list of candidates, to prevent missing GW signals due to the presence of instrumental line-noise closeby.

The search presented in this paper is more robust than but not as sensitive as the hierarchical all-sky search performed by the Einstein@Home distributed computing project on the same data [18], which for example, in the 0.5 Hz-wide band at 152.5 Hz, excluded the presence of signals with a h₀ greater than 7.6 × 10⁻²⁵ at a 90% confidence level. This later run used the Hough transform method to combine the information from coherent searches on a time scale of about a day and it was very computationally intensive. At the same time, the Einstein@Home search, due to its larger coherent baseline, is more sensitive to the fact that the second spin-down is not included in the search. The Hough transform method has also proven to be more robust against transient spectral disturbances than the StackSlide or PowerFlux semi-coherent methods [10].

Other strategies can be applied to perform all-sky multi-interferometer searches using the Hough transform operating on successive SFTs. In this section we estimate the sensitivity of the semi-coherent Hough search for two hypothetical searches to illustrate its capabilities, by either varying the duration of the total amount of data used in the multi-interferometer search, or by lowering the threshold for candidate selection.

In the first case we consider a search over the full S5 data with the same criteria of selecting candidates as presented here, i.e. setting a threshold in the significance for a false alarm level equivalent to one candidate per 0.25 Hz band, but using the full data. This first search would be more sensitive since we increase the number of SFTs to search over. In this case the sensitivity can be estimated from equation (28) and using the desired significance threshold as the 'loudest' event. In this case we should take into account that the number of templates for a two years search is double that for a single year, because of the increase of spin-down values that are resolvable. This corresponds to the dot-dashed line in the significance threshold in figure 6. Different sensitivity confidence levels can also be provided by modifying the false dismissal rate in equation (29) accordingly:

$$\begin{eqnarray*} &&\mathcal{S}_{95\%}= s_{{{\rm th}}}/\sqrt{2} + \textrm {erfc}^{-1}(0.1) , \nonumber \\ &&\mathcal{S}_{90\%}= s_{{{\rm th}}}/\sqrt{2} + \textrm {erfc}^{-1}(0.2) , \nonumber \\ &&\mathcal{S}_{50\%}= s_{{{\rm th}}}/\sqrt{2} + \textrm {erfc}^{-1}(1) . \nonumber \end{eqnarray*}$$

In the second case, we consider only the data from second year of S5 data but lower the threshold in the significance such that the false alarm would be ten candidates per 0.25 Hz band (see dashed line in figure 6).

In figure 14 we show the projected sensitivities for these two searches for different confidence levels. The best sensitivity would be for the search performed on the combined full two years of S5 data. For example in the frequency interval 159.75–160.0 Hz the estimated sensitivity levels are of 8.1 × 10⁻²⁵, 7.9 × 10⁻²⁵ and 7.1 × 10⁻²⁵ at the 95%, 90% and 50% confidence level respectively. In the second case, corresponding to an increase in false alarm rate but with a reduced amount of data, the best sensitivities are of 8.8 × 10⁻²⁵, 8.5 × 10⁻²⁵ and 7.6 × 10⁻²⁵ at the 95%, 90% and 50% confidence level respectively.

Zoom In Zoom Out Reset image size

Notice that although both strategies are, in principle, more sensitive than the search presented in this paper, they would produce many more candidates. These would need eliminating either by demanding coincidence between two searches of comparable sensitivity, or by follow-up using a more sensitive, computationally intensive, search. Coincidence analysis will be explored in future searches owing to the efficiency with which background noise is removed. Follow up studies will always be computationally limited. Therefore the follow up capacity will actually limit the event threshold for candidate selection. Without the inclusion of a coincidence analysis, the event threshold will have to be set higher, therefore compromising the potential sensitivity of the search itself. Moreover, for the first hypothetical case, in order to achieve the projected sensitivity, one would need to perform the search over the entire 67696 SFTs available at once. If one wanted to use the 'look up table' approach to compute the Hough transform over the two years of data, the computational cost would increase by a factor of nine with respect to the one year search presented in this paper, and the memory usage would increase from 0.8 to 7.2 GB, for the same sky-patch size. The memory needs could be reduced by decreasing the sky-patch size, but at additional computational cost. Another consequence of analyzing both years together is that the spin-down step size in the production search would have had to be reduced significantly.

There are a number of areas where further refinements could improve the sensitivity of the Hough search. In particular, one could decrease the grid spacing in parameter space in order to reduce the maximum mismatch allowed, increase the duration of the SFTs to increase the SNR within a single SFT, the development of further veto strategies to increase the overall efficiency of the analysis, as well as the tracking and establishing of appropriate data-cleaning strategies to remove narrow-band disturbances present in the peak-grams [39, 42, 43]. Several of these ideas are being addressed and will be implemented in the 'Frequency Hough all-sky search' using data from the Virgo second and fourth science runs to analyze data between 20 and 128 Hz.

In summary, we have reported the results of an all-sky search for continuous, nearly monochromatic gravitational waves on data from LIGO's fifth science (S5) run, using a new detection pipeline based on the Hough transform. The search covered the frequency range 50–1000 Hz and with the frequency's time derivative in the range −8.9 × 10⁻¹⁰ Hz s⁻¹ to zero. Since no evidence for gravitational waves has been observed, we have derived upper limits on the intrinsic gravitational wave strain amplitude using a standard population-based method. The best upper limits correspond to 1.0 × 10⁻²⁴ for the first year of S5 in the 158–158.25 Hz band, and 8.9 × 10⁻²⁵ for the second year in the 146.5–146.75 Hz band (see figure 10).

This new search pipeline has allowed to process outliers down to significance from 5.10 at 50 Hz to 6.13 at 1000 Hz permitting deeper searchers than in previous Hough all-sky searches [10]. A set of new features have been included into the multi-detector Hough search code to be able to cope with large amounts of data and the memory limitations on the machines. In addition, a χ² veto has been applied for the first time for continuous gravitational wave searches. This veto might be very useful for the analysis of the most recent set of data produced by the LIGO and Virgo interferometers (science runs S6, VSR2 and VSR4) whose data at lower frequencies are characterized by larger contamination of non-Gaussian noise than for S5.

Although the search presented here is not the most sensitive one on the same S5 data, this paper shows the potential of the new pipeline given the advantage of the lower computational cost of the Hough search and its robustness compared to other methods, and suggests further improvements to increase the sensitivity and overall efficiency of the analysis.

The authors gratefully acknowledge the support of the United States National Science Foundation for the construction and operation of the LIGO Laboratory, the Science and Technology Facilities Council of the United Kingdom, the Max-Planck-Society, and the State of Niedersachsen/Germany for support of the construction and operation of the GEO600 detector, and the Italian Istituto Nazionale di Fisica Nucleare and the French Centre National de la Recherche Scientifique for the construction and operation of the Virgo detector. The authors also gratefully acknowledge the support of the research by these agencies and by the Australian Research Council, the International Science Linkages program of the Commonwealth of Australia, the Council of Scientific and Industrial Research of India, the Istituto Nazionale di Fisica Nucleare of Italy, the Spanish Ministerio de Economía y Competitividad, the Conselleria d'Economia Hisenda i Innovació of the Govern de les Illes Balears, the Foundation for Fundamental Research on Matter supported by the Netherlands Organisation for Scientific Research, the Polish Ministry of Science and Higher Education, the FOCUS Programme of Foundation for Polish Science, the Royal Society, the Scottish Funding Council, the Scottish Universities Physics Alliance, The National Aeronautics and Space Administration, the Carnegie Trust, the Leverhulme Trust, the David and Lucile Packard Foundation, the Research Corporation, and the Alfred P Sloan Foundation. This document has been assigned LIGO Laboratory document number LIGO-P1300071.