Implementation of an $\mathcal{F}$-statistic all-sky search for continuous gravitational waves in Virgo VSR1 data

This paper presents results from a wide parameter search for periodic gravitational waves from spinning neutron stars using data from the Virgo detector [1]. The data used in this paper were produced during Virgo’s first science run (VSR1) which started on May 18, 2007 and ended on October 1, 2007. The VSR1 data has never been searched for periodic gravitational waves from isolated neutron stars before. The innovation of the search is the combination of the efficiency of the FFT algorithm together with a nearly optimal grid of templates.

Rotating neutron stars are promising sources of gravitational radiation in the band of ground-based laser interferometric detectors (see [2] for a review). In order for these objects to emit gravitational waves, they must exhibit some asymmetry, which could be due to imperfections of the neutron star crust or the influence of strong internal magnetic fields (see [3] for recent results). Gravitational radiation also may arise due to the r-modes, i.e., rotation-dominated oscillations driven unstable by the gravitational emission (see [4] for discussion of implications of r-modes for GW searches). Neutron star precession is another gravitational wave emission mechanism (see [5] for a recent study). Details of the above mechanisms of gravitational wave emission can be found in [6, 7, 8, 9, 10, 11].

A signal from a rotating neutron star can be modeled independently of the specific mechanism of gravitational wave emission as a nearly periodic wave the frequency of which decreases slowly during the observation time due to the rotational energy loss. The signal registered by an Earth-based detector is both amplitude and phase modulated because of the motion of the Earth with respect to the source.

The gravitational wave signal from such a star is expected to be very weak, and therefore months-long segments of data must be analyzed. The maximum deformation that a neutron star can sustain, measured by the ellipticity parameter $\epsilon$ (Eq. (7)), ranges from $5\times 10^{-6}$ for ordinary matter [7, 12] to $10^{-4}$ for strange quark matter [10, 11]. For unknown neutron stars one needs to search a very large parameter space. As a result, fully coherent, blind searches are computationally prohibitive. To perform a fully coherent search of VSR1 data in real time (i.e., in time of 136 days of duration of the VSR1 run) over the parameter space proposed in this paper would require a $1.4\times 10^{4}$ petaflop computer [13].

A natural way to reduce the computational burden is a hierarchical scheme, where first short segments of data are analyzed coherently, and then results are combined incoherently. This leads to computationally manageable searches at the expense of the signal-to-ratio loss. To perform the hierarchical search presented in this paper in real time a $2.8$ teraflop computer is required. One approach is to use short segments of the order of half an hour long so that the the signal remains within a single Fourier frequency bin in each segment and thus a single Fourier transform suffices to extract the whole power of the gravitational wave signal in each segment. Three schemes were developed for the analysis of Fourier transforms of the short data segments and they were used in the all-sky searches of ground interferometer data: the “stack-slide” [14], the “Hough transform” [15, 14], and the “PowerFlux” methods [14, 16, 17].

Another hierarchical scheme involves using longer segment duration for which signal modulations need to be taken into account. For segments of the order of days long, coherent analysis using the popular $\mathcal{F}$-statistic [18] is computationally demanding, but feasible. For example the hierarchical search of VSR1 data presented in this paper requires around 1900 processor cores for the time of 136 days which was the duration of the VSR1 run. In the hierarchical scheme the first coherent analysis step is followed by a post-processing step where data obtained in the first step are combined incoherently. This scheme was implemented by the distributed volunteer computing project Einstein@Home [19]. The E@H project performed analysis of LIGO S4 and S5 data, leading to results published in three papers. In the first two, [20, 21], the candidate signals from the coherent analysis were searched for coincidences. In the third paper [22] the results from the coherent search were analyzed using a Hough-transform scheme.

The search method used here is similar to the one used in the first two E@H searches: it consists of coherent analysis by means of the $\mathcal{F}$-statistic, followed by a search for coincidences among the candidates obtained in the coherent analysis. There are, however, important differences: in this search we have a fixed threshold for the coherent analysis, resulting in a variable number of candidates from analysis of each of the data segments. Moreover, for different bands we have a variable number of two-day data segments (see Section 3 for details). In E@H searches the number of candidates used in coincidences from each data segment was the same, as was the number of data segments for each band. In addition, the duration of data segments for coherent analysis in the two first E@H searches was 30h (48h in this case).

In this analysis we have implemented algorithms and techniques that considerably improve the efficiency of this search. Most importantly we have used the FFT algorithm to evaluate the $\mathcal{F}$-statistic for two-day data segments. Also we were able to use the FFT algorithm together with a grid of templates that was only 20% denser than the best known grid (i.e., the one with the least number of points).

Given that the data we analyzed (Virgo VSR1) had a higher noise and the duration was shorter than the LIGO S5 data we have achieved a lower sensitivity than in the most recent E@H search [22]. However due to very good efficiency of our code we were able to analyze a much larger parameter space. We have analyzed a large part of the $f-\dot{f}$ plane that was previously unexplored by any of the all-sky searches. For example in this analysis the $f-\dot{f}$ plane searched was $6.5$ larger than that in the full S5 E@H search [22] (see Figure 5). In addition, the data from the Virgo detector is characterized by different spectral artifacts (instrumental and environmental) from these seen in LIGO data. As a results, certain narrow bands excluded from LIGO searches because of highly non-Gaussian noise can be explored in this analysis.

The paper is organized as follows. In Section 2 we describe the data from the Virgo detector VSR1 run, in Section 3 we explain how the data were selected. In Section 4 the response of the detector to a gravitational wave signal from a rotating neutron star is briefly recalled. In Section 5 we introduce the $\mathcal{F}$-statistic. In Section 6 we describe the search method and present the algorithms for an efficient calculation of the $\mathcal{F}$-statistic. In Section 7 we describe the vetoing procedure of the candidates. In Section 8 we present the coincidence algorithm that is used for post-processing of the candidates. Section 9 contains results of the analysis. In Section 10 we determine the sensitivity of this search, and we conclude the paper by Section 11. In A we present a general formula for probability that is used in the estimation of significance of the coincidences.

The VSR1 science run spanned more than four months of data acquisition. This run started on May 18, 2007 at 21:00 UTC and ended on October 1, 2007 at 05:00 UTC. The detector was running close to its design sensitivity with a duty cycle of 81.0% [23]. The data were calibrated, and the time series of the gravitational wave strain $h(t)$ was reconstructed. In the range of frequencies from 10Hz to 10kHz the systematic error in amplitude was 6% and the phase error was 70 mrad below the frequency of 1.9 kHz [24]. A snapshot of the amplitude spectral density of VSR1 data in the Virgo detector band is presented in Figure 1.

The analysis input consists of narrow-band time-domain data segments. In order to obtain these sequences from VSR1 data, we have used the software described in [25], and extracted the segments from the Short Fourier Transform Database (SFDB). The time domain data sequences were extracted with a certain sampling time $\Delta t$, time duration ${T_{\mathrm{obs}}}$ and offset frequency $f_{\rm off}$. Thus each time domain sequence has the band $[f_{\rm off},f_{\rm off}+B]$, where $B=1/2\times 1/\Delta t$. We choose the segment duration to be exactly two sidereal days and the sampling time equal to $1/2$ s, i.e., the bandwidth $B=1$ Hz. As a result each narrow band time segment contains $N=344\,656$ data points. We have considered 67 two-day time frames for the analysis. The starting time of the analysis (the time of the first sample of the data in the first time frame) was May 19, 2007 at 00:00 UTC. The data in each time frame $d=1,\dots,67$ were divided into narrow band sequences of 1Hz band each. The bands overlap by $2^{-5}$ Hz resulting in 929 bands numbered from 0 to 928. The relation between the band number $b$ and the offset frequency $f_{\rm off}$ is thus given by

Consequently, we have 62 243 narrow band data segments. From this set we selected good data using the following conditions. Let $N_{0}$ be the number of zeros in a given data segment (data point set to zero means that at that time there was no science data). Let $[f_{\rm off},f_{\rm off}+B]$ be the band of a given data segment. Let $S_{min}$ and $S_{max}$ be the minimum and the maximum of the amplitude spectral density in the interval $[f_{\rm off}+0.05B,f_{\rm off}+B-0.05B]$. Spectral density was estimated by dividing the data of a given segment into short stretches and averaging spectra of all the short stretches. We consider the data segment as good data and use it in this analysis, if the following two criteria are met:

1. 1.

   $N_{0}/N<1/4$,

2. 2.

   $S_{max}/S_{min}\leq 1.1$.

20 419 data segments met the above two criteria. Figure 2 shows a time-frequency distribution of the good data segments.

The good data appeared in 50 out of 67 two-day time frames and in 785 out of 929 bands. The distribution of good data segments in time frames and in bands is given in Figure 3.

In Figure 4 we plot a snapshot spectral density of VSR1 data presented in Figure 1 and an estimate of the spectral density of the data that was used in the analysis. We estimate the spectral density in each band by a harmonic mean of the spectral densities of each of the two-day time segments chosen for the analysis in the band. As the data spectrum in each band is approximately white, we estimate the spectral density in each time segment by $2\sigma^{2}\Delta t$, where $\sigma^{2}$ is the variance of the data in a segment and $\Delta t$ is the $1/2$s sampling time. We plot the band from 100Hz to 1 kHz which we have chosen for this search.

The dimensionless noise-free response $h$ of a gravitational-wave detector to a weak plane gravitational wave in the long wavelength approximation, i.e., when the characteristic size of the detector is much smaller than the reduced wavelength $\lambda/(2\pi)$ of the wave, can be written as the linear combination of the two independent wave polarizations $h_{+}$ and $h_{\times}$,

where $F_{+}$ and $F_{\times}$ are the detector’s beam-pattern functions,

The beam-patterns $F_{+}$ and $F_{\times}$ are linear combinations of $\sin 2\psi$ and $\cos 2\psi$, where $\psi$ is the polarization angle of the wave. The functions $a(t)$ and $b(t)$ are the amplitude modulation functions, that depend on the location and orientation of the detector on the Earth and on the position of the gravitational-wave source in the sky, described in the equatorial coordinate system by the right ascension $\alpha$ and the declination $\delta$ angles. They are periodic functions of time with the period of one and two sidereal days. The analytic form of the functions $a(t)$ and $b(t)$ for the case of interferometric detectors is given by Eqs. (12) and (13) of [18]. For a rotating nonaxisymmetric neutron star, the wave polarization functions are of the form

where $h_{0+}$ and $h_{0\times}$ are constant amplitudes of the two polarizations and $\phi(t)+\phi_{0}$ is the phase of the wave, $\phi_{0}$ being the initial phase of the waveform. The amplitudes $h_{0+}$ and $h_{0\times}$ depend on the physical mechanism responsible for the gravitational radiation, e.g., if a neutron star is a triaxial ellipsoid rotating around a principal axis with frequency $f$, then these amplitudes are

where $\iota$ is the angle between the star’s angular momentum vector and the direction from the star to the Earth, and the amplitude $h_{0}$ is given by

Here $I$ is the star’s moment of inertia with respect to the rotation axis, $r$ is the distance to the star, and $\epsilon$ is the star’s ellipticity defined by $\epsilon=|I_{1}-I_{2}|/I$, where $I_{1}$ and $I_{2}$ are moments of inertia with respect to the principal axes orthogonal to the rotation axis. We assume that the gravitational waveform given by Eqs. (2)–(5) is almost monochromatic around some angular frequency $\omega_{0}$, which we define as instantaneous angular frequency evaluated at the solar system barycenter (SSB) at $t=0$, and we assume that the frequency evolution is accurately described by one spindown parameter $\omega_{1}$. Then the phase $\phi(t)$ is given by

where, neglecting the relativistic effects, ${\bf r}_{d}(t)$ is the vector that joins the SSB with the detector, and ${\bf n}$ is the unit vector pointing from SSB to the source. In equatorial coordinates $(\delta$, $\alpha)$ we have ${\bf n}=(\cos\delta\cos\alpha,\cos\delta\sin\alpha,\sin\delta)$.

A method to search for gravitational wave signals from a rotating neutron star in a detector data $x(t),t=1,...,N$ uses the $\mathcal{F}$-statistic, described in [18]. The $\mathcal{F}$-statistic is obtained by maximizing the likelihood function with respect to the four unknown parameters - $h_{0}$, $\phi_{0}$, $\iota$, and $\psi$. This leaves a function of only the remaining four parameters - $\omega_{0}$, $\omega_{1}$, $\delta$, and $\alpha$. Thus the dimension of the parameter space that we need to search decreases from 8 to 4. In this analysis we shall use an observation time ${T_{\mathrm{obs}}}$ equal to the integer multiple of sidereal days. Since the bandwidth of the signal over our coherent observation time of two days is very small, we can assume that over this band the spectral density of the noise is white (constant). Under these assumptions the $\mathcal{F}$-statistic is given by [26, 27]

where $\sigma^{2}$ is the variance of the data, and

The search consists of two parts; the first part is a coherent search of two-day data segments, where we search a 4-parameter space defined by angular frequency $\omega_{0}$, angular frequency derivative $\omega_{1}$, declination $\delta$, and right ascension $\alpha$. The search is performed on a 4-dimensional grid in the parameter space described in Section 6.2. We set a fixed threshold of 20 for the $\mathcal{F}$-statistic for each data segment. This corresponds to a threshold of 6 for the signal-to-noise ratio. All the threshold crossings are recorded together with corresponding 4 parameters of the grid point and the signal-to-noise ratio $\rho$. The signal-to-noise is calculated from the value of the $\mathcal{F}$-statistic at the threshold crossing as

In this way for each narrow band segment we obtain a set of candidates. The candidates are then subject to the vetoing procedure described in Section 7. The second part of the search is the post-processing stage involving search for coincidences among the candidates. The coincidence procedure is described in Section 8.

We apply three vetoing criteria to the candidates obtained in the coherent part of the search - line width veto, stationary line veto and polar caps veto. Data from the detector always contain some periodic interferences (lines) that are detector artifacts. An important part of our vetoing procedure was to identify the lines in the data. We have therefore performed a Fourier search with frequency resolution of 1/(2 days) $\simeq 5.8\times 10^{-6}$ Hz for periodic signals of each of the two-day data segments. We compared the frequencies of the significant periodic signals identified by our analysis with the line frequencies obtained by the Virgo LineMonitor and we found that all the lines from the LineMonitor were detected by our Fourier search.

In order to find coincidences among the candidates, we applied a method similar to the one used in the first two E@H searches [20, 21]. For each band we have searched for coincidences among candidates in different time frames. We are able to search for coincidences only in those bands where there were two or more time frames with data selected for the analysis. If we search for a real gravitational wave signal, we must take into account frequency evolution due to spindown of the rotating neutron star. Thus the first step in the coincidence analysis was to transform all frequencies $\omega_{0}(t_{l})$ of the candidates to a common fiducial reference time $t_{f}$. We have chosen the fiducial time to be the time of the first sample of the latest time frame that we analyzed i.e., the 67th time frame. We have

where $t_{l}$ is the time of the first sample of the $l$th time frame. The next step was to divide the parameter space into cells. This construction of the coincidence cell was different from that in the E@H analysis. To construct the cells in the parameter space we have used the reduced Fisher matrix $\tilde{\Gamma}$ for the linear signal model defined by Eqs. (17) and (18). The reduced Fisher matrix is the projected Fisher matrix on the 4-dimensional space spanned by parameters $\bf{\kappa}=(\omega_{0},\omega_{1},\alpha_{1},\alpha_{2})$. We define the cell in the parameter space by the condition:

Because the ephemeris of the detector is different in each of the time frames, the reduced Fisher matrix is different in each time frame. To have a common coincidence grid we have chosen the grid defined by the latest frame i.e., the frame no. 67 as the coincidence grid. After the transformation of the candidate frequencies to a reference time and construction of the coincidence grid, the coincidence algorithm for each of the bands proceeded in the following steps:

1. 1.

   Transform angles $\alpha$ and $\delta$ to $\alpha_{1}$ and $\alpha_{2}$ coordinates (see Eq.(19))

2. 2.

   Transform candidate parameters to $x_{l}$ coordinates defined by

   $$x_{l}=\sum_{k=1}^{4}\kappa_{k}V_{kl}\sqrt{e_{l}},$$

   (28)

   where $V_{kl},k=1,..,4$ are eigenvectors of the matrix $\tilde{\Gamma}$, and $e_{l}$ are its eigenvalues. In these coordinates the Fisher matrix is proportional to the unit matrix.

3. 3.

   Coordinates $x_{l}$ are rounded to the nearest integer. In this way we sort candidates efficiently into adjacent 4-dimensional hypercubes. If there are more than one candidates from a given data segment in a hypercube we select the candidate that has the highest SNR. We do sorting for each time frame in the band. If there is more than one candidate in a given hypercube we register a coincidence.

4. 4.

   We shift cubes by 1/2 of their size in all possible $2^{4}$ directions, and for each shift we search for coincidences.

This last step of the algorithm takes into account cases for which the candidate events are located on opposite sides of cell borders, edges, and corners and consequently coincidences that could not be found just by packing candidates into adjacent cells.

The most significant coincidence in each band is the one which has the highest multiplicity. For each most significant coincidence we have calculated the false alarm probability i.e., the probability that such a coincidence can occur purely by chance. The false alarm probability is calculated using the formula explained in A and given by Eq. (35). This general formula applies to a variable number of candidates in various time slots and also takes into account the $2^{4}$ shifts of the cells in the parameter space.

In this analysis we have searched coherently 20 419 two-day time segments of data narrowbanded to $1$ Hz. In the coherent part of the search described in Section 6 we have used $9.10\times 10^{16}$ templates which is the number of $\mathcal{F}$-statistic values computed. This resulted in 20 419 candidate files containing $4.21\times 10^{10}$ candidates. The candidates were subject to the vetoing using the three veto criteria: line veto, polar caps veto, and stationary line veto described in Section 7. As a result of vetoing around 24% of the candidates were discarded leaving $3.19\times 10^{10}$ candidates. Nearly all candidates were vetoed by the line veto, whereas 0.20% were vetoed by the stationary line criterion and only $3.2\times 10^{-2}$% by the polar caps veto. In Figure 7 we present an example of the candidate distribution obtained from the coherent search of one narrow band data segment and after the vetoing procedure.

In the next step we have searched for the significant coincidences among the candidates. We have searched for coincidences in all the frequency bands where there were two or more data segments analyzed. In Figure 8 we have plotted the highest coincidence multiplicity for each of the bands. The highest multiplicity was 6, and it occurred in 10 bands. The multiplicity tends to grow with the frequency, because the size of the parameter space grows as $f^{3}$.

For each band we have calculated the false alarm probability corresponding to the most significant coincidence using Eq. (35). The most significant coincidence occurred in band no. 401, corresponding to the frequency range of $\sim[488.4,489.4]$ Hz. It was a coincidence of multiplicity = 5, and its false alarm probability was 14.5%. Note that a coincidence with the highest multiplicity is not the most significant. This is because the significance depends on the number of time frames with candidates in a given band and also on the number of candidates in the time frames. By adopting a criterion used by E@H searches that the background coincidences correspond to false alarm probability of 0.1% or greater, we conclude that we have found no significant coincidence and thus no viable gravitational wave candidate. Considering the significance of the coincidences, we could adopt even a 10% false alarm probability as a background. Consequently, we proceed to the final stage of our analysis - estimation of sensitivity of the search.

The sensitivity of the search is taken to be the amplitude $h_{0}$ of the gravitational wave signal that can be confidently detected. To estimate the sensitivity we use a procedure developed in [20]. We determine the sensitivity of the search in each of the 785 frequency bands that we have searched. To determine the sensitivity, we perform Monte Carlo simulations in which, for a given amplitude $h_{0}$, we randomly select the other seven parameters of the signal: $\omega_{0},\omega_{1},\alpha,\delta,\phi_{0},\iota,$ and $\psi$. We choose frequency and spindown parameters uniformly over their range, and source positions uniformly over the sky. We choose angles $\phi_{0}$ and $\psi$ uniformly over the interval $[0,2\pi]$ and we choose $\cos\iota$ uniformly over the interval $[-1,1]$. For each band we add the signal to all the data segments chosen for the analysis in that band. Then we process the data through our pipeline. First, we perform a coherent $\mathcal{F}$-statistic search of each of the data segments where the signal was added, and store all the candidates above our $\mathcal{F}$-statistic threshold of 20. In this coherent analysis, to make the computation manageable, we search over only parameter space consisting of $\pm$ 2 grid points around the nearest grid point where the signal was added. Then we apply our vetoing procedure to the candidates obtained as explained in Section 7. Finally, we perform coincidence analysis of the candidates that survive vetoing which is described in Section 8. We define a detectable signal if it is coincident in more than 70% of the time frames in a given band. This condition is similar to the condition used in the two E@H searches, where a coincidence method was used [20, 21]. For bands with only one frame available the coherent search over one $2$-day data segment was performed. In this case the injected signal is declared detected if its signal-to-noise ratio obtained in the coherent search is larger than the signal-to-noise ratio of the loudest signal in that data segment without an injection. For each band we inject signals with 5 different amplitude values, and perform 100 randomized injections for each amplitude. For each amplitude we calculate how many signals were detected, and by interpolation we determine the amplitude corresponding to 90% of signals detected. This amplitude was defined as the 90% confidence sensitivity. Sometimes even for the highest amplitude we have not reached the $90\%$ detection probability. In this case we performed injections for higher amplitudes until the desired level of detectability was achieved. In Figure 9, as an example, we present an estimation of the sensitivity in the band $b=369$ corresponding to the frequency range of $\sim[457.47,458.44]$ Hz. The errors in the sensitivity estimates originate from calibration errors in the amplitude and errors due to a finite number of Monte Carlo injections. We use 100 injections; hence from a binomial statistic, one $\sigma$ is equivalent to 3% fluctuation. Thus the estimated amplitude sensitivity corresponds to confidence in the range from 87% to 93%. To estimate how this uncertainty in confidence translates into uncertainty in the amplitude we have performed an additional set of injections for a range of amplitudes close to the estimated sensitivity and from the slope of the confidence vs. amplitude we determined the uncertainty in the amplitude. To increase the accuracy of the error estimate, we have performed 1000 injections for each amplitude. The uncertainty in the amplitude was not more than 5%. The calibration errors in VSR1 data are 6% (see Section 2). Adding these two types of errors in quadrature results in the total error in sensitivity estimate to be around 7%.

The sensitivity of this search obtained through Monte Carlo simulations for the whole band searched is presented in Figure 10.

We see from Figure 10 that the sensitivity essentially reflects the instrumental noise curve given in Figure 4. We have made a fit of the sensitivity $h^{90\%}_{0}$ to the one-sided spectral density $S_{h}(f)$ of detector noise described by the following relation:

We find that the prefactor $R_{\cal D}$ is in the range from 15.6 to 22.4 and it depends on the frequency band and the number of the data segments in the band.

The sensitivity of this search was 50% to 2 times better, depending on the bandwidth, than that of the LIGO S4 search [20] and comparable to the sensitivities obtained in the early LIGO S5 data [21], but $2$ to $5$ times worse than the upper limits in the E@H full LIGO S5 data search [22]. This was due to the lower noise and longer observation time for the LIGO S5 data w.r.t the Virgo VSR1 data. However, for the first time in an all-sky search we have estimated the sensitivity in the frequency band from 400Hz to 1kHz and the frequency spindown range from $-7.2\,(f_{0}/400\,{\rm Hz})\times 10^{-9}\,$ Hz/s to $-6.0\times 10^{-9}\,$ Hz/s, which is a previously unexplored region in the parameter space.

The next step is to test the search method described in this paper in the Mock Data Challenge (MDC) designed by the LIGO and Virgo projects to validate and compare pipelines that are proposed to be used in the analysis of the forthcoming data from the advanced detectors. It is also planned to further test the pipeline presented here with other data sets collected by the LIGO and Virgo detectors.