All-sky search of NAUTILUS data
A search for periodic gravitational-wave signals from isolated neutron stars in the NAUTILUS detector data is presented. We have analyzed half a year of data over the frequency band Hz, the spindown range Hz/s and over the entire sky. We have divided the data into 2 day stretches and we have analyzed each stretch coherently using matched filtering. We have imposed a low threshold for the optimal detection statistic to obtain a set of candidates that are further examined for coincidences among various data stretches. For some candidates we have also investigated the change of the signal-to-noise ratio when we increase the observation time from two to four days. Our analysis has not revealed any gravitational-wave signals. Therefore we have imposed upper limits on the dimensionless gravitational-wave amplitude over the parameter space that we have searched. Depending on frequency, our upper limit ranges from to . We have attempted a statistical verification of the hypotheses leading to our conclusions. We estimate that our upper limit is accurate to within %.
We present results of the search of the NAUTILUS (a resonant bar detector ) data for periodic gravitational-wave (GW) signals. We search the band Hz and the spindown range Hz/s over the entire sky. We analyze NAUTILUS data collected in the year 2001. We divide the data into stretches of 2 sidereal days. Each stretch of data is analyzed coherently using matched filtering in the form of the -statistic [2, 3]. We have analyzed slightly more than half a year of data.
Previous analysis of bar detector data for periodic GW signals were the search of the galactic center and the globular cluster 47 Tucanae with the ALEGRO detector , the search of the galactic center using the EXPLORER detector data , and an all-sky search using the EXPLORER data [6, 7]. LIGO (Laser Gravitational Wave Observatory) data was searched for known pulsars [8, 9, 10], over the entire sky  using the coherent method and over all the sky using incoherent methods [12, 13]. Currently LIGO data are analyzed over the entire sky by the Einstein@Home project .
In section 2 we present the data analysis methods used in our search. In section 3 we outline our search procedure. In section 4 we discuss the analysis of the candidates. This analysis consists of two parts: the first part is the search for coincidences between the candidates each obtained from a different stretch of data and the second part is an investigation of the increase of the signal-to-noise ratio of candidates when we increase the observation time from two to four days. In section 5 we impose upper limits on amplitudes of the gravitational waves in the parameter space that we have searched.
2 Data analysis methods
In order to search for gravitational waves from long lived periodic sources we have used the maximum likelihood (ML) method. For the case of Gaussian noise the ML method consists of linearly filtering the data with a template matched to the signal that we are searching for. The main complication of the matched filtering is that the signal depends on several unknown parameters. This requires evaluation of the likelihood function over a large parameter space. In order to minimize the computation time we use several data analysis tools. Firstly we find the maximum likelihood estimators of some parameters (4 in the case of a GW signal from a rotating neutron star that we are searching for) in a closed analytic form, thereby reducing the dimensionality of the parameter space that we have to search. The likelihood function over the reduced parameter space is called the -statistic and it is derived in . Secondly we analyze data of length equal to an integer multiple of a sidereal day. This leads to a considerable simplification of the -statistic and consequently reduces the number of numerical operations to evaluate it. The -statistic for observation time equal to an integer number of sidereal days is given in section III of . Thirdly we use optimal numerical algorithms, in particular the Fast Fourier Transform (FFT) in order to calculate the -statistic efficiently. Fourthly we minimize the number of -statistic calculations over the parameter space by solving a covering problem for this space [16, 17]. Let us explain the latter two data analysis tools in more detail. The response of a bar detector to a gravitational-wave signal from a spinning neutron star is summarized in section 2.1 of .
Fast Fourier Transform. Estimates have shown [21, 20] that for the bandwidth and the spindown range that we search we need to take into account in our templates only one spindown parameter in order to match the signal. Consequently the phase modulation function of the waveform is given by
where is angular frequency and is the spindown parameter, is the constant unit vector in the direction of the star in the Solar System Barycenter (SSB) reference frame (it depends on the right ascension and the declination of the source), and is the vector joining the SSB with the detector and is the speed of light. The detection statistic involves two integrals of the form
where is the data stream, is the amplitude modulation function that depends on and . The above integral is not a Fourier transform because the frequency in the phase multiplies the term which is a non linear function of time. In order to convert the integral into a Fourier transform we introduce the following interpolation procedure. The phase [Equation (1)] can be written as
The functions and do not depend on the angular frequency . We can write the integral (2) as
We see that the integral (6) can be interpreted as a Fourier transform (and computed efficiently with an FFT), if . In order to convert equation (6) to a Fourier transform we introduce a new time variable , so called barycentric time [18, 2],
Thus in order to compute the integral (6), we first multiply the data by the function for each set of the parameters and then resample the resulting function according to equation (7). At the end we perform the FFT.
The covering problem. The covering problem is to find the minimum number of templates in the parameter space , so that the fractional loss in signal to ratio is not less than ( is the minimal match parameter introduced by Owen ). In order to solve the covering problem we introduce a useful approximate model of the gravitational-wave signal from a rotating neutron star. The model relies on (i) neglecting all spindowns in the phase modulation due to motion of the detector with respect to the SSB; and (ii) discarding the component of the vector (connecting the SSB and the detector) which is perpendicular to the ecliptic plane. These approximations lead to the following phase model of the signal:
where and are new constant parameters,
where is the obliquity of the ecliptic and where and are known functions of time,
is the -component of the vector joining the center of Earth and the SSB, and is the -component of the vector joining the center of Earth and the detector. and are the corresponding -components. We also neglect the slowly varying modulation of the signal’s amplitude, so finally we approximate the whole signal by
where and are the constant amplitude and initial phase, respectively. The above signal model is called linear because it has the property that its phase given by equation (9) is a linear function of the parameters. We have shown  that the above model is a good approximation to the accurate response of the detector to the GW signal in the sense that the Fisher matrix for the linear model reproduces well the Fisher matrix for the accurate model. Thus whenever a Fisher matrix is needed we can use the Fisher matrix for the linear model as an approximation to the Fisher matrix for the accurate model. The great advantage of the linear model is that components of its Fisher matrix are constant, independent of the values of the parameters. In order to solve the covering problem for the parameter space we use the Fisher matrix as a metric on the parameter space. Because the components of the Fisher matrix are constant the grid is uniform what greatly simplifies its construction. In our search, as a grid we use the hypercubic lattice . However we have an additional constraint. In order to apply the FFT algorithm the nodes of the grid must coincide with the Fourier frequencies. We have constructed a suitable grid by performing rotations and dilatations of the original hypercubic lattice. We construct the grid in the parameters and then transform it to parameters for which the -statistic is calculated.
The linear parametrization has one more application. We use it in order to calculate the threshold for the -statistic corresponding to a certain false alarm probability. Namely, using the linear parametrization we divide the parameter space into cells as explained in [2, 20]. All the cells are exactly the same and their number is easily calculated using the Fisher matrix (see section IIIB of ). The false alarm probability is the probability that exceeds threshold in one or more cells and is given by
where is the false alarm probability for a single cell.
3 Search procedure
We have searched the data collected by the NAUTILUS detector in the year 2001. The bandwidth of Hz, where the detector is most sensitive, has been analyzed. We have divided the data into stretches which span two sidereal days. We have assumed a minimum pulsar spindown age equal to yrs and so we have searched the negative frequency time derivatives in the range of Hz/s. For this and two days of the observation time it is sufficient to include only one spindown in the phase [21, 20]. Each two day sequence was analyzed coherently using the -statistic. We have used the constrained hypercubic grid as explained in the previous section. For the grid construction we have assumed the minimal match parameter . Using this minimal-match value our grid consists of around grid points ( frequency bins, spindowns, sky positions). The threshold on 2 corresponding to false alarm probability has been calculated using Equation (15) and is around 72. In order to compensate the loss of signal-to-noise ratio (SNR) due to the discreteness of the grid, imperfect templates and numerical approximation in evaluation of the -statistic (resampling procedure) we have adopted two lower thresholds on 2 equal to and . We have registered parameters of all templates which crossed the threshold of . For threshold crossings of we have performed a verification procedure. The verification procedure consisted of calculating the -statistic for the template parameters of the candidate using a four day stretch of data involving the original two day stretch. For a true gravitational-wave signal by this procedure one would expect an increase of signal-to-noise ratio by . In total, we have analyzed 93 two day data stretches. In Figure 1 we have presented the two-sided amplitude spectrum of the NAUTILUS detector data that we have analyzed. The spectrum was obtained in the following way. We have estimated the power spectrum density in each of the 93 two day data sequences and then we have taken the square root of the average of the 93 power spectra. We see that the best sensitivity is around Hz. Moreover we have obtained the rms error of our power spectrum estimate by calculating the variance from the estimates of the spectra of in each of the 93 data segments. The relative 1 error in the amplitude power spectrum is around 18%.
During the search we have obtained candidates above the -threshold of and above the -threshold of .
4 Analysis of the candidates
4.1 Signal-to-noise ratio of the candidates
In Figure 2 we have plotted a histogram of the frequencies of all the candidates above the -threshold of 50. The histogram shows an excess of candidates in the frequency band of Hz. This excess is a result of the presence of a periodic interference in the data that appears as a series of harmonics in the bandwidth of the detector. One of the harmonics is located in the subband Hz. The effect of the harmonic is visible in our estimate of the spectrum (Figure 1) and appears as a bump in the band Hz.
As a first step in the candidate analysis we have calculated the increase in signal-to-noise ratio when we increase the observation time from two to four days. This has been done for all the candidates above the threshold . Figure 3 shows the highest increase in SNR for candidates when going from a two day data stretch to the four day one. The maximum is calculated for each of the 93 data stretches analyzed.
We see that typically the highest gain in the signal-to-noise is . This should be compared with the theoretical gain of of SNR when we increase the observation time by a factor of 2. The periodic interference present in the data to which we attribute these maximum SNR increases does not gives a higher increase of the SNR because its frequency changes erratically over the observation time of days and it cannot reproduce the Doppler shift of a real GW signal modulated by detector motion with respect to the SSB. Assuming that the two day sequence is independent of the four day sequence we could perform the -test that consists of calculation of the ratio of the -statistic for 4 days observation time and the -statistic for 2 days observation time. Taking as the null hypothesis for the test that data is only Gaussian noise the -statistic has the central distribution with 4-degrees of freedom and the ratio has Fisher-Snedeckor distribution . The typical highest value of for a given data segment is around The probability of crossing the threshold is around %. This would give a high confidence that data is noise only. Unfortunately this is only a crude approximation because the two day sequence is contained in the four day one and the assumption of independence of the two -statistic is not fulfilled.
4.2 Coincidences among the candidates from different data stretches
Candidates from different data stretches are considered coincident if they cluster closely together in the four-dimensional parameter space . We employ the clustering method described in , which uses a grid of “coincidence cells”. This method will reliably detect strong signals which would produce candidates with closely-matched parameters in many of the different data stretches.
In a first step, the frequency value of each candidate above the threshold of is shifted to the same fiducial time: the GPS start time of the earliest () stretch, . Defining to be the time span of two sidereal days, the frequencies of the candidates are shifted to via
where is the starting time of the ’th data stretch, given by .
To find coincidences, a grid of cells is constructed such that the cells are rectangular in the coordinates . The dimensions of the cells are adapted to the parameter space search. Thus, the cells are constructed to be as small as possible to reduce the probability of coincidences due to false alarms. However, since each of the 93 different data stretches uses a slightly different parameter space grid, the coincidence cells must be chosen to be large enough that the candidates from a source (which would appear at slightly different points in parameter space in each of the 93 data stretches) would still lie in the same coincidence cell. As a conservative choice we use cell sizes in of , in of , and an isotropic cell grid in the sky with equatorial spacing of . Each candidate event is assigned to a particular cell. In cases where two or more candidate events from the same data stretch fall into the same cell, only the candidate having the largest value of is retained in the cell. Then the number of candidate events per cell coming from distinct data stretches is counted.
From the 93 different data stretches, this coincidence method found that we get candidates which appear consistently in no more than 4 data stretches uniformly over the search bandwidth, where there are no instrumental interferences. This is the background of the number of coincidences. We would like to test the null hypothesis that the coincidences are result of the noise only. Let us assume that the parameter space is divided into independent coincidence cells, the candidate events are independent and the probability for a candidate event to fall into any given coincidence cell is . Thus probability that a given coincidence cell is populated with one or more candidate event is given by
where is the number of candidate events per data segment. The probability that any given coincidence cell contains candidate events from or more distinct data segments is given by a binomial distribution
Finally the probability that there is or more coincidences in one or more of the cells is
The expected number of cells with or more coincidences is given by
In our case the number of cells is given by , the number of data segments is , and the number of candidates per data segment is . From Equation (19) we find that the probability of finding or more coincident candidates is almost one. Thus for the background coincidences we can accept the null hypothesis that they are from noise only with a high confidence. Over the bandwidth Hz we find an excess of coincidences with the maximum of 8 coincidences. By Equation (19), the false alarm probability associated with with 8 or more coincidences is of the order of and thus they cannot be attributed to noise. We consider these coincidences to be due to the periodic interference present in the data.
5 Upper limits
Our verification procedure consisting of coincidences among the candidates from distinct data segments and an analysis of the increase of signal-to-noise ratio presented in section 4 did not produce convincing evidence of a gravitational-wave signal. We therefore proceeded to estimate the upper limits for the amplitudes of the gravitational-wave signals in the parameter space that we have searched. Detection of a signal is signified by a large value of the -statistic that is unlikely to arise from the noise-only distribution. If instead the value of is consistent with pure noise with high probability we can place an upper limit on the strength of the signal. One way of doing this is to take the loudest event obtained in the search and solve the equation
for signal-to-noise ratio , where is the detection probability, is the value of the -statistic corresponding to the loudest event, and is a chosen confidence. Then is the desired upper limit with confidence . We can also obtain an upper limit with confidence for several independent searches from the equation
where is the threshold corresponding to the loudest event in ’th search and is the number of searches. Here is the probability that a signal with signal-to-noise ratio crosses the threshold in at least one of the independent searches. To calculate we assume that the data have a Gaussian distribution and consequently the probability of detection has a non-central distribution with 4 degrees of freedom and the noncentrality parameter equal to . We have investigated this assumption by obtaining histograms of the -statistic values of the candidates and comparing them to the central distribution with 4 degrees of freedom. The result is shown in Figure 4. There is an overall qualitative agreement of candidates distributions with the theoretical one. However, the candidates distributions do not pass a goodness-of-fit test for a distribution at the significance level of 5%.
In order to translate our upper limit on the SNR into the upper limit on the gravitational-wave amplitude, we use the Equation (93) of  for signal-to-noise ratio of a GW signal from a spinning neutron star averaged over the source position and orientation. Thus and are related by the following formula:
where is one-sided spectral density at frequency . We have used Equations (22) and (23) to obtain upper limits in Hz bands over the bandwidth Hz that we have searched. The upper limit results are presented in Figure 5. Assumed Gaussian noise, we have chosen the confidence and we denote the upper limits by .
Our best upper limit is equal to at a frequency of 922.55 Hz. Using our 1 rms error of the amplitude power spectrum estimate we reckon that our upper limit has likewise an error of %.
We would like to thank Maria Alessandra Papa for making available to us large computing resources of the Albert-Einstein-Institut enabling to complete the data analysis for this paper in a reasonable time. We would also like to thank the Interdisciplinary Center for Mathematical and Computational Modelling of Warsaw University for computing time. A Królak would like to acknowledge the financial support of Istituto Nazionale di Fisica Nucleare (INFN), the INFN Section in University of Rome “La Sapienza” and the Max-Planck-Institut für Gravitationsphysik (Albert-Einstein-Institut). H Pletsch would like to acknowledge the support of the IMPRS on Gravitational Wave Astronomy.
-  Astone P et al(ROG Collaboration) 1997 Astroparticle Physics 7 231-243
-  Jaranowski P, Królak A and Schutz B F 1998 Phys. Rev. D 58 063001
-  Jaranowski P and Królak A 2005 Living Rev. Relativity 8 3
-  Mauceli et al2000 arXiv:gr-qc/0007023v1
-  Astone P et al2002 Phys. Rev. D 65 022001
-  Astone P et al2003 Class. Quantum Grav. 20 S665–S676
-  Astone P et al2005 Class. Quantum Grav. 22 S1243–S1254
-  Abbott B et al(LIGO Scientific Collaboration) 2004 Phys. Rev. D 69 082004
-  Abbott B et al(LIGO Scientific Collaboration), Kramer M and Lyne A 2005 Phys. Rev. Lett. 94 181103
-  Abbott B et al(LIGO Scientific Collaboration), Kramer M and Lyne A 2007 Phys. Rev. D 76 042001
-  Abbott B et al(LIGO Scientific Collaboration) 2007 Phys. Rev. D 76 082001
-  Abbott B et al(LIGO Scientific Collaboration) 2005 Phys. Rev. D 72 102004
-  Abbott B et al(LIGO Scientific Collaboration) 2008 Phys. Rev. D 77 022001
-  LIGO Scientific Collaboration 2008 The Einstein@Home search for periodic gravitational waves in LIGO S4 data Preprint arXiv:0804.1747v1 [gr-qc]
-  Astone P, Borkowski K M, Jaranowski P and Królak A 2002 Phys. Rev. D 65 042003
-  Conway J H and Sloane N J A 1993 Sphere Packings, Lattices and Groups (New York: Springer-Verlag)
-  Prix R 2007 Class. Quantum Grav.24 S481–S490
-  Schutz B F 1991 The Detection of Gravitational Waves ed D G Blair (Cambridge: Cambridge University Press) pp 406–452
-  Jaranowski P and Królak A 1999 Phys. Rev. D 59 063003
-  Jaranowski P and Królak A 2000 Phys. Rev. D 61 062001
-  Brady P R, Creighton T, Cutler C and Schutz B F 1998 Phys. Rev. D 57 2101
-  Owen B 1996 Phys. Rev. D 53 6749