Improved Hough search for gravitational wave pulsars 6 Alicia M. Sintes†‡, and Badri Krishnan‡ 0 † Departament de F´ısica, Universitat deles Illes Balears, Cra. Valldemossa Km. 7.5, E-07122 0 Palma de Mallorca, Spain 2 ‡ Max-Planck-Institut fu¨r Gravitationsphysik, Albert Einstein Institut,Am Mu¨hlenberg 1, n D-14476 Golm, Germany a J E-mail: [email protected], [email protected] 9 1 Abstract. We describe an improved version of the Hough transform search for continuous gravitational waves from isolated neutron stars assuming the input to be short segments 1 of Fourier transformed data. The method presented here takes into account possible non- v stationarities of the detector noise and the amplitude modulation due to the motion of the 1 detector. Thesetwoeffectsaretakenintoaccountforthefirststageonly,i.e. thepeakselection, 8 tocreate thetime-frequencymap of ourdata,while theHoughtransform itself isperformed in 0 1 the standard way. 0 6 0 / c q 1. Introduction - TheHoughtransformisapatternrecognition algorithmwhichwasoriginallyinventedtoanalyze r g bubblechamberpicturesfromCERN[1]. ItwaslaterpatentedbyIBM[2],andithasfoundmany : v applications in the analysis of digital images [3]. A detailed discussion of the Hough transform i as applied to the search for continuous gravitational waves can be found in [4]. An example of X such a search is [5] in which this semi-coherent technique is used to perform an all-sky search r a for isolated spinning neutron stars using two months of data collected in early 2003 from the second science run of the LIGO detectors [6, 7]. The main results of [5] are all-sky upper limits on a set of narrow frequency bands within the range 200-400Hz and including one spin-down parameter. The best upper limit on the gravitational wave strain amplitude that we obtained in this frequency range is 4.43×10−23. Several searches have been completed or are underway using the same data. Examples of such searches are [8] which targets known radio pulsars at twice thepulsar frequency, and [9] that performsan all sky search in a widefrequency band, but using only about 10 hours of data and applying coherent (matched filtering) techniques. Such techniques are very computationally intensive for large parameter space searches. The ultimate goal for wide parameter space searches for continuous signals over large data sets is to employ hierarchical schemes [10, 11, 12, 13, 14, 15] which alternate coherent and semi- coherent techniques. The Hough transform search would then be used to select candidates in parameter space to be followed up. It is crucial that the first stage of the search is as efficient as possible, since it determines the final sensitivity of the full pipeline. Inthis paperwepresent animproved implementation of theHough transformwith respectto what was described in [4] when this is applied to a set of Fourier transformed data (without any demodulations). Thisnewimplementationtakesintoaccountnon-stationaritiesofthenoisefloor and also the amplitude modulation due to the antenna pattern that were previously ignored. The plan of the paper is as following: Section 2 described the waveforms we are looking for. Section 3 describes the basics of the standard Hough transform and in Section 4 how this could be improved in the presence of non-stationarities. 2. The signal at the detector The output of the detector can be represented by x(t)= h(t)+n(t), (1) where n(t) is the noise that affects the observation at a detector time t, and h(t) is the gravitational wave signal received at the detector: h(t) = F (n,ψ)h (t)+F (n,ψ)h (t). (2) + + × × Here h and h are the two independent polarizations of the strain amplitude at the detector + × and F and F are the antenna beam pattern functions, that depend on the direction n to the + × star and also on the polarization angle ψ. Due to the motion of the Earth, the F depend +,× implicitly on time. In case of a laser interferometer detector with perpendicular arms, the expressions for F and F are given by [16]: + × F = a(t)cos2ψ+b(t)sin2ψ, F =b(t)cos2ψ−a(t)sin2ψ, (3) + × where the functions a(t) and b(t) are independent of ψ. We shall make the standard assumption that over a certain time-scale, the noise n(t) is stationary, although in realistic cases this hypothesis is likely to be violated at some level. Within this approximation, the Fourier components n˜(f) of the noise are statistically described by: 1 E[n˜(f)n˜∗(f′)] = δ(f −f′)S (f), (4) n 2 where E[] denotes the expectation value with respect to an ensemble of noise realization, the ∗ superscript denotes complex conjugate, S (f) is the one sided noise power spectral density, and n tildes denote Fourier transforms. The form of the gravitational wave emitted by an isolated pulsar depends on the physical mechanism which may cause the neutron star to emit periodic gravitational waves. The main possibilities considered in the literature are (i) non-axisymmetric distortions of the solid part of the star, (ii) unstable r-modes in the fluid, and (iii) free precession (or ‘wobble’) (see [17] and references therein). If we assumethe emission mechanism is (i), then thegravitational waves are emittedatafrequencywhichistwicetherotationalratef ofthepulsar. Underthisassumption, r the waveforms for the two polarizations h are given by [16]: +,× 1+cos2ι h = A cosΦ(t) =h cosΦ(t), h = A sinΦ(t) = h cosιsinΦ(t), (5) + + 0 × × 0 2 whereι is the angle between the pulsar’sspin axis andthe direction of propagation of the waves, and h is the amplitude: 0 16π2GI ǫf2 h = zz r . (6) 0 c4 d Here G is Newton’s gravitational constant, c the speed of light, I is the z-z component of the zz star’s moment of inertia with the z-axis being its spin axis, ǫ := (I −I )/I is the equatorial xx yy zz ellipticity of the star, and d is the distance of the star from Earth. The phase Φ(t) takes its simplest form in the Solar System Barycenter (SSB) frame where it can be expanded in a Taylor series: f Φ(t)= Φ +2π (n) (T −T )n+1. (7) 0 0 (n+1)! nX≥0 Here T is time in the SSB frame and T is a fiducial start time. The phase Φ , frequency f 0 0 (0) and the spin-down parameters f , n > 0, are defined at this fiducial start time. Neglecting (n) relativistic effects which do not affect us significantly in this case, the relation between the time of arrival T of the wave in the SSB frame and in the detector frame t is n·r T = t+ , (8) c where n is the direction to the pulsar and r is the detector position in the SSB frame. The instantaneous frequency f(t) of the wave as observed by the detector is given, to a very good approximation, by the familiar non-relativistic Doppler formula: v(t)·n f(t)−fˆ(t) =fˆ(t) (9) c where v(t) is the velocity of the detector at the detector time t, fˆ(t) is the instantaneous signal frequency at time t: f fˆ(t) = f + (n)(T −T )n. (10) (0) 0 n! nX>0 Equations (9) and (10) describe the time-frequency pattern produced by a signal, and this is the pattern that the Hough transform is used to look for. 3. Standard Hough transform The Hough transform can be used to find the pattern produced by the Doppler shift (9) and the spin-down (10) of a gravitational wave signal in the time-frequency plane of our data. The parameters which determine this pattern are ({f },n). The parameter space is covered by a (n) discrete cubic grid and the result of the Hough transform is an histogram for each point of this grid. By “standard” implementation of the Hough transform, we refer to the implementation usedin [5]wherethestartingpointareN shortsegments of Fourier transformeddata(which are called Short Fourier Transforms, or SFTs) with T being the time baseline of the SFTs. From coh this set of SFTs, calculating the periodograms (the square modulus of the Fourier transform) and selecting frequency bins (peaks) above a certain threshold ρ , we obtain a time-frequency th map of our data. Assuming a fixed threshold, the optimal value turns out to be ρ = 1.6, th corresponding to a peak selection probability of q = e−ρth = 0.2 in the absence of a signal [4]. For each selected bin in the SFTs, we find which points in parameter space are consistent with it, according to Eq. (9), and the number count in all such points is increased by unity. This is repeated for all the selected bins in all the SFTs to obtain the final histogram. There are several criteria to select the frequency bins. One of them is to select by setting a threshold ρ on the normalized power ρ defined as th k 2|x˜ |2 k ρ = , (11) k T S (f ) coh n k where x˜ is the discrete Fourier transform of the data, the frequency index k corresponds to a k physical frequency of f = k/T , and S (f ) is the single sided power spectral density of the k coh n k detector noise. The kth frequency bin is selected if ρ ≥ ρ , and rejected otherwise. In this k th way, each SFT is replaced by a collection of zeros and ones called a peak-gram. The Hough transform is used to calculate the number count n at each parameter space point starting from this collection of peak-grams, i.e. each point in parameter space corresponds to a pattern in the time-frequency plane of our data, and the number count n is the sum of the ones and zeros of the different peak-grams along this curve. Let p(n) be the probability distribution of n in the absence of a signal, and p(n|h) the distribution in the presence of a signal h(t). It is clear that 0 ≤ n≤ N, where N is the number of SFTs, and it can be shown that for stationary Gaussian noise, p(n) is a binomial distribution with mean Nq where q = e−ρth is the probability that any frequency bin is selected: N p(n)= qn(1−q)N−n. (12) (cid:18) n (cid:19) In the presence of a signal, the distribution p(n|h) is ideally also a binomial but with a slightly larger mean Nη where, for weak signals, η is given by ρ η = q 1+ thλ +O(λ2) . (13) 2 k k n o λ is the signal to noise ratio within a single SFT, and for the case when there is no mismatch k between the signal and the template: 4|h˜(f )|2 k λ = (14) k T S (f ) coh n k with h˜(f) being the Fourier transform of the signal h(t) (see [4] for details) Candidates in parameter space are selected by setting a threshold n on the number count. th For the case of large N, and Gaussian stationary noise, the relation between the the false alarm rate α and the threshold n is independent of the signal strength and given by: th n = Nq+ 2Nq(1−q)erfc−1(2α) (15) th p where, erfc−1 is the inverse of the complementary error function. On the average, the weakest signal which will cross the above thresholds at a false alarm rate α and false dismissal β is given by S1/2 S h = 5.34 n , S = erfc−1(2α)+erfc−1(2β). (16) 0 N1/4rT coh Equation (16) gives the smallest signal which can be detected by the search. 4. Improved peak selection for the Hough transform The approximation that the distribution p(n|h) in the presence of a signal is binomial can break down for several reasons: the random mismatch between the signal and the template used to calculate the number count, non-stationarity of the noise, and the amplitude modulation of the signal which causes λ to vary from one SFT to another and for different sky locations and k pulsar orientations. The result of these effects is to “smear” out the binomial distribution in the presence of a signal as it is illustrated in Figure 14 in [5]. To account for these non-stationarities a possible Adaptive Hough Transform is discussed in [18] in which Palomba et al suggest not to increase the number-count by unity, but by a real number depending on the detector response function and the noise level. In this paper we propose an alternative method that consists in setting an adaptive peak selection threshold that varies for the different SFTs and sky locations, and then performing the Standard Hough Transform on the resulting peak-grams. We expect this method should have a somewhat similar gain in sensitivity as the method presented in [18], but with a better performance since all the implementation tricks described in [4] and [5] Sec VI.C can still be applied. However, we emphasize that a proper statistical justification of this method is still unclear and is being investigated. 4.1. Dealing with non-stationarity Let us consider the quantity ρ λ /2, related to the second summand in right hand side of k k Eq. (13), which, in the discrete case, can be rewritten as: ρ λ |x˜ |2 |h˜(f )|2 |x˜ |2 hE[|n˜ |2]i |h˜(f )|2 k k k k k k N k = = . (17) 2 E[|n˜ |2] E[|n˜ |2] E[|n˜ |2] E[|n˜ |2] hE[|n˜ |2]i k k k k k N Here E[|n˜ |2] = T S (f )/2 denotes the noise contribution in a single SFT, which is assumed k coh n k to be stationary in a time scale T although it can vary from one SFT to another. In practice coh the expected noise level, in each data stretch, is estimated from the data itself by means of the running median [19] and it is frequency dependent. hE[|n˜ |2]i denotes the average noise floor k N of the N different SFTs. WetypicallychooseT = 30min. WhilealargervalueofT wouldleadtobettersensitivity, coh coh this time baseline cannot be made arbitrarily large because of the frequency drift caused by the Doppler effect and the spin-down; we would like the signal power of a putative signal to be concentrated in less than half the frequency resolution 1/T . It is shown in [4] that at 300Hz, coh we could ideally choose T up to ∼ 60min. On the other hand, we should be able to find a coh significantnumberofsuchdatastretchesduringwhichtheinterferometersareinlock,thenoiseis stationary, and the data are labeled satisfactory according to certain data quality requirements. In order to keep ρ λ /2 constant across the different N data segments, Eq. (17) suggests th k that, ignoring the amplitude modulation, a more natural threshold for the peak selection should be placed not on ρ but on k |x˜ |2 hE[|n˜ |2]i k k N ρˆ = . (18) k E[|n˜ |2] E[|n˜ |2] k k The probability of selecting a peak in the absence of a signal in a single SFT would then be: E[|n˜ |2] k qˆ = exp −ρ (19) k (cid:20) th hE[|n˜ |2]i (cid:21) k N thus selecting fewer peaks in the most noisy periods. The presence of very noisy periods, if those are not vetoed for the analysis, could in principle disturb not just some SFTs but also dominate the average hE[|n˜ |2]i , and this could degrade k N the sensitivity. Therefore, in practice, hE[|n˜ |2]i should be estimated by a more sophisticated k N method,e.g. basedon themean ofsome sortedpoints between differentquartiles of thedifferent SFTs, or the harmonic mean etc. 4.2. Amplitude modulation corrections In Eq. (14) the amplitude of the signal |h˜(f )|2 is proportional to F2A2 +F2A2, again taking k + + × × a and b to be approximately constant. Averaging over the polarization angle ψ 1 2π h···i := dψ(···), (20) ψ 2π Z 0 we have hF2i = hF2i = (a2 + b2)/2. Therefore h|h˜(f )|2i ∝ a2 + b2 and this holds true + ψ × ψ k ψ independently of the polarization of the waveform, and also |h˜(f )|2 ∝ a2 + b2 for circular k polarization (cosι= 1) independently of ψ. The functions a and b vary in time due to the motion of the detector, although they can be considered constant for the typical time baseline T of the SFTs. This suggests that the peak coh selection criteria should be based on setting a threshold on: |x˜ |2 hE[|n˜ |2]i a2+b2 k k N ρ¯ = , (21) k E[|n˜ |2] E[|n˜ |2] ha2+b2i k k N andselecting thosefrequencybins(peaks)whichsatisfyρ¯ > ρ . Thepeakselection probability k th in the absence of a signal would be: E[|n˜ |2] ha2+b2i k N q¯ = exp −ρ . (22) k (cid:20) th hE[|n˜ |2]i a2+b2 (cid:21) k N In other words, the threshold on |x˜ |2/E[|n˜ |2] would be frequency dependent and also different k k for different SFTs and sky locations. Since the functions a and b vary for different sky locations and we want to take advantage of the implementation method described in [4], for a wide area or an all-sky search, the celestial sphere should then be divided into several smaller regions. For each of these small sky-patches, the functions a and b should be evaluated at the center of the sky-patch and assumed constant for the entire sky-patch. The peak selection procedure needs to be repeated for each individual sky-patch we want to analyze. The smaller the patch the more the gain in sensitivity and also the smaller the distortions induced by the tiling of the celestial as discussed in [5]. The details of the statistical properties of this search will be presented elsewhere. In particular, as mentioned earlier, a proper statistical justification of this method in terms of, say, a Neyman-Pearson criteria, the regime where the gain in sensitivity is most significant, and a large scale Monte-Carlo study are currently being investigated. Acknowledgments We would like to thank the LIGO Scientific Collaboration for many useful discussions. We also acknowledge the support of the Max-Planck Society and the Spanish Ministerio de Educacio´n y Ciencia Research Project REF: FPA-2004-03666. AMS is grateful to the Albert Einstein Institute for hospitality where this work was initiated. Bibliography [1] Hough P V C 1959 in Proceedings of the International Conference on High Energy Accelerators and Instrumentation, edited by L Kowarski, CERN, p.554 [2] Hough P V C 1962 U.S. Patent 3,069,654 [3] Illingworth J and KittlerJ 1988 Computer Vision, Graphics, and Image Processing 44 87-116 [4] Krishnan B, Sintes A M, Papa M A, SchutzB F, Frasca S and Palomba C 2004 Phys. Rev. D 70 082001 [5] AbbottB et al (TheLIGO ScientificCollaboration) 2005 Phys. Rev. D 72 102004 [6] AbramoviciA et al 1992 Science 256 325 [7] Barish B and Weiss R 1999 Phys. Today 52 (10) 44 [8] AbbottB et al (TheLIGO ScientificCollaboration) 2005 Phys. Rev. Lett. 94 181103 [9] AbbottB et al (TheLIGO ScientificCollaboration) 2005 LIGO-P050008-00-Z [10] Brady P R and Creighton T 2000 Phys. Rev. D 61 082001 [11] PapaMA,SchutzBFandSintesAM2001,inGravitational waves: Achallenge totheoretical astrophysics, ICTP Lecture Notes Series, Vol. III,edited byV. Ferrari, J.C. Miller, L. Rezzolla, Italy,p. 431 [12] Cutler, C Gholami I and Krishnan B 2005 Phys. Rev. D 72 042004 [13] Frasca S 2000 Int. J. Mod. Phys. D 9 369 [14] Frasca S, AstoneP and Palomba C 2005 Class. Quantum Grav. 22 S1013 [15] Frasca S and Palomba C 2004 Class. Quantum Grav. 21 S1645 [16] Jaranowski P, Kr´olak A and SchutzB F1998 Phys. Rev. D 58 063001 [17] AbbottB et al (TheLIGO ScientificCollaboration) 2004 Phys. Rev. D 69 082004 [18] Palomba C, AstoneP and Frasca S 2005 Class. Quantum Grav. 22 S1255 [19] Mohanty S D 2002 Class. Quantum Grav. 19 1513