\mathcal{F}-statistic all-sky search for continuous gravitational waves in Virgo data

Implementation of an F-statistic all-sky search for continuous gravitational waves in Virgo VSR1 data

Abstract

We present an implementation of the -statistic to carry out the first search in data from the Virgo laser interferometric gravitational wave detector for periodic gravitational waves from a priori unknown, isolated rotating neutron stars. We searched a frequency range from 100 Hz to 1 kHz and the frequency dependent spindown range from  Hz/s to zero. A large part of this frequency - spindown space was unexplored by any of the all-sky searches published so far. Our method consisted of a coherent search over two-day periods using the -statistic, followed by a search for coincidences among the candidates from the two-day segments. We have introduced a number of novel techniques and algorithms that allow the use of the Fast Fourier Transform (FFT) algorithm in the coherent part of the search resulting in a fifty-fold speed-up in computation of the -statistic with respect to the algorithm used in the other pipelines. No significant gravitational wave signal was found. The sensitivity of the search was estimated by injecting signals into the data. In the most sensitive parts of the detector band more than 90% of signals would have been detected with dimensionless gravitational-wave amplitude greater than .

pacs:
04.80.Nn, 95.55.Ym, 97.60.Gb, 07.05.Kf

1 Introduction

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 (Eq. (7)), ranges from for ordinary matter [7, 12] to 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 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 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 -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 -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 -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 plane that was previously unexplored by any of the all-sky searches. For example in this analysis the plane searched was 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 -statistic. In Section 6 we describe the search method and present the algorithms for an efficient calculation of the -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.

2 Data from the Virgo’s first science run

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 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.

3 Data selection

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 , time duration and offset frequency . Thus each time domain sequence has the band , where . We choose the segment duration to be exactly two sidereal days and the sampling time equal to  s, i.e., the bandwidth  Hz. As a result each narrow band time segment contains 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 were divided into narrow band sequences of 1Hz band each. The bands overlap by  Hz resulting in 929 bands numbered from 0 to 928. The relation between the band number and the offset frequency is thus given by

 foff[Hz]=100+(1−2−5)b. (1)

Consequently, we have 62 243 narrow band data segments. From this set we selected good data using the following conditions. Let 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 be the band of a given data segment. Let and be the minimum and the maximum of the amplitude spectral density in the interval . 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. ,

2. .

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 , where is the variance of the data in a segment and is the s sampling time. We plot the band from 100Hz to 1 kHz which we have chosen for this search.

4 The response of the detector

The dimensionless noise-free response 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 of the wave, can be written as the linear combination of the two independent wave polarizations and ,

 h(t)=F+(t)h+(t)+F×(t)h×(t), (2)

where and are the detector’s beam-pattern functions,

 F+(t)=a(t)cos2ψ+b(t)sin2ψ, (3) F×(t)=b(t)cos2ψ−a(t)sin2ψ. (4)

The beam-patterns and are linear combinations of and , where is the polarization angle of the wave. The functions and 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 and the declination angles. They are periodic functions of time with the period of one and two sidereal days. The analytic form of the functions and 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

 h+(t)=h0+cos(ϕ(t)+ϕ0),h×(t)=h0×sin(ϕ(t)+ϕ0), (5)

where and are constant amplitudes of the two polarizations and is the phase of the wave, being the initial phase of the waveform. The amplitudes and 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 , then these amplitudes are

 h0+=12h0(1+cos2ι),h0×=h0cosι, (6)

where is the angle between the star’s angular momentum vector and the direction from the star to the Earth, and the amplitude is given by

 h0=16π2Gc4ϵIf2r. (7)

Here is the star’s moment of inertia with respect to the rotation axis, is the distance to the star, and is the star’s ellipticity defined by , where and 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 , which we define as instantaneous angular frequency evaluated at the solar system barycenter (SSB) at , and we assume that the frequency evolution is accurately described by one spindown parameter . Then the phase is given by

 ϕ(t)=ω0t+ω1t2+n⋅rd(t)c(ω0+2ω1t), (8)

where, neglecting the relativistic effects, is the vector that joins the SSB with the detector, and is the unit vector pointing from SSB to the source. In equatorial coordinates , we have .

5 The F - statistic

A method to search for gravitational wave signals from a rotating neutron star in a detector data uses the -statistic, described in [18]. The -statistic is obtained by maximizing the likelihood function with respect to the four unknown parameters - , , , and . This leaves a function of only the remaining four parameters - , , , and . 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 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 -statistic is given by [26, 27]

 F≈2σ2(|Fa|2⟨a2⟩+|Fb|2⟨b2⟩), (9)

where is the variance of the data, and

 Fa:=N∑t=1x(t)a(t)exp[−iϕ(t)], (10) Fb:=N∑t=1x(t)b(t)exp[−iϕ(t)].
 ⟨a2⟩=N∑t=1a(t)2,⟨b2⟩=N∑t=1b(t)2. (11)

6 Description of the search

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 , angular frequency derivative , declination , and right ascension . 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 -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 . The signal-to-noise is calculated from the value of the -statistic at the threshold crossing as

 ρ=√2(F−2). (12)

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.

6.1 Choice of the parameter space

We have searched the frequency band from 100 Hz to 1 kHz over the entire sky. We have followed [13] to constrain the maximum value of the parameter for a given frequency by , where is the minimum spindown age1. We have chosen yr for the whole frequency band searched. Also, in this search we have considered only the negative values for the parameter , thus assuming that the rotating neutron star is spinning down. This gives the frequency-dependent range of the spindown parameter where :

 |f1|≤1.6×10−9f100Hz1000yrτmin[Hzs−1] (13)

We have considered only one frequency derivative. Estimates taking into account parameter correlations (see [13] Figure 6 and Eqs. (6.2) - (6.6)) show that even for the minimum spindown age of yr and for two days coherent observation time that we consider here, it is sufficient to include just one spindown parameter.

In Figure 5 we have compared the parameter space searched in this analysis in the plane with that of other recently published all-sky searches: Einstein@Home early S5 search [21], Einstein@Home full S5 [22], PowerFlux early S5 [28], PowerFlux full S5 [17].

6.2 Efficient calculation of the F-statistic on the grid in the parameter space

Calculation of the -statistic (Eq. (9)) involves two sums given by Eqs. (10). By introducing a new time variable called the barycentric time ([29, 18, 27])

 tb:=t+n⋅rd(t)c. (14)

we can write these sums as discrete Fourier transforms in the following way

 Fa=N∑tb=1x(t(tb))a(t(tb))exp[−iϕs(t(tb))]exp[−iω0tb], (15) Fb=N∑tb=1x(t(tb))b(t(tb))exp[−iϕs(t(tb))]exp[−iω0tb],

where

 ϕs(t)=ω1t2+2n⋅rd(t)cω1t. (16)

Written in this form, the two sums can be evaluated using the Fast Fourier Transform (FFT) algorithm thus speeding up their computation dramatically. The time transformation described by equation (14) is called resampling. In addition to the use of the FFT algorithm we apply an interpolation of the FFT using the interbinning procedure (see [27] Section VB). This results in the -statistic sampled twice as fine with respect to the standard FFT. This procedure is much faster than the interpolation of the FFT obtained by padding the data with zeroes and calculating a FFT that is twice as long. With the approximations described above for each value of the parameters , , and , we calculate the -statistic efficiently for all the frequency bins in the data segment of bandwidth 1Hz.

In order to search the 4-dimensional parameter space, we need to construct a 4-dimensional grid. To minimize the computational cost we construct a grid that has the smallest number of points for a certain assumed minimal match [30]. This problem is equivalent to the covering problem [31, 32] and it has the optimal solution in 4-dimensions in the case of a lattice i.e., a uniformly spaced grid. In order that our parameter space is a lattice, the signals’ reduced Fisher matrix must have components that are independent of the values of the parameters. This is not the case for the signal given by Eqs. (2) - (8); it can be realized however for an approximate signal called the linear model described in Section IIIB of [27]. The linear model consists of neglecting amplitude modulation of the signal and discarding the component of the vector joining the detector and the solar system barycenter that is perpendicular to the ecliptic. This approximation is justified because the amplitude modulation is very slow compared to the phase modulation and the discarded component in the phase is small compared to the others. As a result the linear model signal has a constant amplitude , and one can find parameters such that the phases are linear functions of them. We explicitly have (see Section IIIB of [27] for details):

 hlin(t)=A0cos[ϕlin(t)+ϕ0], (17)

where

 ϕlin(t)=ω0t+ω1t2+α1μ1(t)+α2μ2(t). (18)

The parameters and are defined by

 α1:=ω0(sinαcosδcosε+sinδsinε), (19) α2:=ω0cosαcosδ, (20)

where is the obliquity of the ecliptic, and and are known functions of the detector ephemeris.

In order that the grid is compatible with application of the FFT, its points should be constrained to coincide with Fourier frequencies at which the FFT is calculated. Moreover, we observed that a numerically accurate implementation of the interpolation to the barycentric time (see Eq. (14)) is so computationally demanding that it may offset the advantage of the FFT. Therefore we introduced another constraint in the grid such that the resampling is needed only once per sky position for all the spindown values. Construction of the constrained grid is described in detail in Section IV of [27]. In this search we have chosen the value of the minimal match . The workflow of the coherent part of the search procedure is presented in Figure 6.

7 Vetoing procedure

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)  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.

7.1 Line width veto

We veto all the candidates with frequency around every known line frequency according to the following criterion

 |f−fl|<Δfmax, (21)

where the width is estimated as