Data analysis for continuous gravitational wave signals

1
Data analysis for continuous gravitational wave
signals

• Alicia M. Sintes
• Universitat de les Illes Balears, Spain

Villa Mondragone International School of
Gravitation and Cosmology Frascati, September 9th
2004
Photo credit NASA/CXC/SAO
2
Talk overview

• Gravitational waves from pulsars
• Sensitivity to pulsars
• Target searches
• Review of the LSC methods results
• Frequency domain method
• Time domain method
• Hardware injection of fake pulsars
• Pulsars in binary systems
• The problem of blind surveys
• Incoherent searches
• The Radon Hough transform
• Summary and future outlook

3
Gravitational waves from pulsars

• Pulsars (spinning neutron stars) are known to
  exist!
• Emit gravitational waves if they are
  non-axisymmetric

4
GWs from pulsars The signal

• The GW signal from a neutron star
• Nearly-monochromatic continuous signal
• spin precession at frot
• excited oscillatory modes such as the r-mode at
  4/3 frot
• non-axisymmetric distortion of crystalline
  structure, at 2frot

• (Signal-to-noise)2

5
Signal received from a pulsar

• A gravitational wave signal we detect from a
  pulsar will be
• Frequency modulated by relative motion of
  detector and source
• Amplitude modulated by the motion of the antenna
  pattern of the detector

6
Signal received from an isolated pulsar

• The detected strain has the form

strain antenna patterns of the detector to plus
and cross polarization, bounded between -1 and 1.
They depend on the orientation of the detector
and source and on the polarization of the waves.
the two independent wave polarizations
the phase of the received signal depends on the
initial phase and on the frequency evolution of
the signal. This depends on the spin-down
parameters and on the Doppler modulation, thus on
the frequency of the signal and on the
instantaneous relative velocity between source
and detector
T(t) is the time of arrival of a signal at the
solar system barycenter (SSB).Accurate timing
routines (2?s) are needed to convert between GPS
and SSB time. Maximum phase mismatch 10-2
radians for a few months
7
Signal model isolated non-precessing neutron
star
Limit our search to gravitational waves from an
isolated tri-axial neutron star emitted at twice
its rotational frequency (for the examples
presented here, only)

• h0 - amplitude of the gravitational wave signal
• ? - angle between the pulsar spin axis and line
  of sight

- equatorial ellipticity
8
Sensitivity to Pulsars
9
Target known pulsars

• For target searches only one search template (or
  a reduced parameter space) is required
• e.g. search for known radio pulsars with
  frequencies (2frot) in detector band . There are
  38 known isolated radio pulsars fGW gt 50 Hz
• Known parameters position, frequency and
  spin-down (or approximately)
• Unknown parameters amplitude, orientation,
  polarization and phase
• Timing information provided from radio
  observations
• In the event of a glitch we would need to add an
  extra parameter for jump in GW phase
• Coherent search methods can be used
• i.e. those that take into account amplitude and
  phase information

10
Crab pulsar

• Young pulsar (60 Hz) with large spin-down and
  timing noise. Frequency residuals will cause
  large deviations in phase if not taken into
  account.
• E.g., use Jodrell Bank monthly crab ephemeris to
  recalculate spin down parameters. Phase
  correction can be applied using interpolation
  between monthly ephemeris.

11
Review of the LSC methods and results

• S1 run Aug 23 - Sep 9 2002 (400 hours).LIGO,
  plus GEO and TAMA
• Setting upper limits on the strength of periodic
  gravitational waves using the first science data
  from the GEO600 and LIGO detectors Phys. Rev.
  D69, 082004 (2004), gr-qc/0308050,

GEO LLO-4K LHO-4K LHO-2K 3x Coinc.
Duty cycle 98 42 58 73 24
12
Directed Search in S1
NO DETECTION EXPECTED at present sensitivities
Detectable amplitudes with a 1 false alarm rate
and 10 false dismissal rate Upper limits on from
spin-down measurements of known radio pulsars
Crab Pulsar
Predicted signal for rotating neutron star with
equatorial ellipticity e d I/I 10-3 , 10-4 ,
10-5 _at_ 8.5 kpc.
PSR J19392134 1283.86 Hz P 1.0511 10-19
s/s D 3.6 kpc
13
Two search methods

• Frequency domain
• Conceived as a module in a hierarchical search
• Best suited for large parameter space
  searches(when signal characteristics are
  uncertain)
• Straightforward implementation of standard
  matched filtering technique (maximum likelihood
  detection method)
• Cross-correlation of the signal with the
  template and inverse weights with the noise
• Frequentist approach used to cast upper limits.

• Time domain
• process signal to remove frequency variations due
  to Earths motion around Sun and spindown
• Best suited to target known objects, even if
  phase evolution is complicated
• Efficiently handless missing data
• Upper limits interpretation Bayesian approach

14
Frequency domain search

• The input data are a set of SFTs of the time
  domain data, with a time baseline such that the
  instantaneous frequency of a putative signal does
  not move by more than half a frequency bin.The
  original data being calibrated, high-pass
  filtered windowed.
• Data are studied in a narrow frequency band.
  Sh(f) is estimated from each SFT.
• Dirichlet Kernel used to combine data from
  different SFTs (efficiently implements matched
  filtering, or other alternative methods).
• Detection statistic used is described in
  Jaranoski, Krolak, Schutz, Phys. Rev.
  D58(1998)063001
• The F detection statistic provides the maximum
  value of the likelihood ratio with respect to the
  unknown parameters,
  , given the data and the template
  parameters that are known

15
Frequency domain method

• The outcome of a target search is a number F
  that represents the optimal detection statistic
  for this search.
• 2F is a random variable For Gaussian stationary
  noise, follows a c2 distribution with 4 degrees
  of freedom with a non-centrality parameter
  l?(hh). Fixing ?, ? and ?0 , for every h0, we
  can obtain a pdf curve p(2Fh0)
• The frequentist approach says the data will
  contain a signal with amplitude ? h0 , with
  confidence C, if in repeated experiments, some
  fraction of trials C would yield a value of the
  detection statistics ? F
• Use signal injection Monte Carlos to measure
  Probability Distribution Function (PDF) of F

16
Measured PDFs for the F statistic with fake
injected worst-case signals at nearby frequencies
h0 1.9E-21
h0 2.7E-22
Note hundreds of thousands of injections were
needed to get such nice clean statistics!
95
95
2F
2F
2F 1.5 Chance probability 83
2F 3.6 Chance probability 46
h0 5.4E-22
h0 4.0E-22
95
95
2F
2F
2F 6.0 chance probability 20
2F 3.4 chance probability 49
17
Time domain target search

• Method developed to handle known complex phase
  evolution. Computationally cheap.
• Time-domain data are successively heterodyned to
  reduce the sample rate and take account of pulsar
  slowdown and Doppler shift,
• Coarse stage (fixed frequency) 16384 ? 4
  samples/sec
• Fine stage (Doppler spin-down correction) ? 1
  samples/min ? Bk
• Low-pass filter these data in each step. The data
  is down-sampled via averaging, yielding one value
  Bk of the complex time series, every 60 seconds
• Noise level is estimated from the variance of the
  data over each minute to account for
  non-stationarity. ? ?k
• Standard Bayesian parameter fitting problem,
  using time-domain model for signal -- a function
  of the unknown source parameters h0 ,?, ? and ?0

18
Time domain Bayesian approach

• We take a Bayesian approach, and determine the
  joint posterior distribution of the probability
  of our unknown parameters, using uniform priors
  on h0 ,cos ?, ? and ?0 over their accessible
  values, i.e.
• The likelihood ? exp(-?2 /2), where
• To get the posterior PDF for h0, marginalizing
  with respect to the nuisance parameters cos ?, ?
  and ?0 given the data Bk

19
Upper limit definition detection
The 95 upper credible limit is set by the value
h95 satisfying Such an upper limit can be
defined even when signal is present A detection
would appear as a maximum significantly offset
from zero
20
Posterior PDFs for CW time domain analyses
Simulated injection at 2.2 x10-21
p
shaded area 95 of total area
p
21
S1 Results

• No evidence of CW emission from PSR J19392134.
• Summary of 95 upper limits for ho

IFO Frequentist FDS Bayesian TDS GEO
(1.9?0.1) x 10-21 (2.2?0.1) x 10-21 LLO
(2.7?0.3) x 10-22 (1.4?0.1) x 10-22
LHO-2K (5.4?0.6) x 10-22 (3.3?0.3) x
10-22 LHO-4K (4.0?0.5) x 10-22
(2.4?0.2) x 10-22
Previous results for this pulsar ho lt 10-20
(Glasgow, Hough et al., 1983), ho lt 1.5 x
10-17 (Caltech, Hereld, 1983).
22
S2 Upper LimitsFeb 14 Apr 14, 2003
95 upper limits

• Performed joint coherent analysis for 28
  pulsars using data from all IFOs
  
  
  
  
  
  
• Most stringent UL is for pulsar J1910-5959D
  (221 Hz) where 95 confident that h0 lt
  1.7x10-24
• 95 upper limit for Crab pulsar ( 60 Hz) is h0
  lt 4.7 x 10-23
• 95 upper limit for J19392134 ( 1284 Hz) is h0
  lt 1.3 x 10-23

23
Equatorial Ellipticity

• Results on h0 can be interpreted as upper limit
  on equatorial ellipticity
• Ellipticity scales with the difference in radii
  along x and y axes

• Distance r to pulsar is known, Izz is assumed to
  be typical, 1045 g cm2

24
S2 hardware injections

• Performed end-to-end validation of analysis
  pipeline by injecting simultaneous fake
  continuous-wave signals into interferometers
• Two simulated pulsars were injected in the LIGO
  interferometers for a period of 12 hours during
  S2
• All the parameters of the injected signals were
  successfully inferred from the data

25
Preliminary results for P1

• Parameters of P1

P1 Constant Intrinsic Frequency Sky position
0.3766960246 latitude (radians) 5.1471621319
longitude (radians) Signal parameters are defined
at SSB GPS time 733967667.026112310 which
corresponds to a wavefront passing LHO at GPS
time 733967713.000000000 LLO at GPS time
733967713.007730720 In the SSB the signal is
defined by f 1279.123456789012 Hz fdot 0 phi
0 psi 0 iota p/2 h0 2.0 x 10-21
26
Preliminary results for P2

• Parameters for P2

P2 Spinning Down Sky position 1.23456789012345
latitude (radians) 2.345678901234567890
longitude (radians) Signal parameters are defined
at SSB GPS time SSB 733967751.522490380, which
corresponds to a wavefront passing LHO at GPS
time 733967713.000000000 LLO at GPS time
733967713.001640320 In the SSB at that moment the
signal is defined by f1288.901234567890123 fdot
-10-8 phase2 pi (f dt1/2 fdot
dt2...) phi 0 psi 0 iota p/2 h0 2.0 x
10-21
27
GW pulsars in binary systems

• Physical scenarios
• Accretion induced temperature asymmetry
  (Bildsten, 1998 Ushomirsky, Cutler, Bildsten,
  2000 Wagoner, 1984)
• R-modes (Andersson et al, 1999 Wagoner, 2002)
• LMXB frequencies are clustered (could be detected
  by advanced LIGO).
• We need to take into account the additional
  Doppler effect produced by the source motion
• 3 parameters for circular orbit
• 5 parameters for eccentric orbit
• possible relativistic corrections
• Frequency unknown a priori

28
Coincidence Pipeline
H1 SFTs
L1 SFTs
Compute F statistic over bank of filters and
frequency range Store results above threshold
Compute F statistic over bank of filters and
frequency range Store results above threshold
Consistent events in parameter space?
no
yes
FL1 gt Fchi2 FH1 gt Fchi2
FL1 lt Fchi2 FH1 lt Fchi2
FL1 or FH1 lt Fchi2
Pass chi2 cut ?
Pass chi2 cut?
no
Candidates
Rejected events
29
All-Sky and targeted surveys for unknown pulsars

• It is necessary to search for every signal
  template distinguishable in parameter space.
  Number of parameter points required for a
  coherent T107s search
• Brady et al., Phys.Rev.D57 (1998)2101
• Number of templates grows dramatically with the
  integration time. To search this many parameter
  space coherently, with the optimum sensitivity
  that can be achieved by matched filtering, is
  computationally prohibitive.
• gtIt is necessary to explore alternative search
  strategies

Class f (Hz) t (Yrs) Ns Directed All-sky
Slow-old lt200 gt103 1 3.7x106 1.1x1010
Fast-old lt103 gt103 1 1.2x108 1.3x1016
Slow-young lt200 gt40 3 8.5x1012 1.7x1018
Fast-young lt103 gt40 3 1.4x1015 8x1021
30
Alternative search strategies

• The idea is to perform a search over the total
  observation time using an incoherent
  (sub-optimal) method
• We propose to search for evidence of a signal
  whose frequency is changing over time in
  precisely the pattern expected for some one of
  the parameter sets
• The methods used are
• Radon transform
• Hough transform
• Power-flux method
• Phase information is lost between data segments

31
The Radon transform

• Break up data into N segments
• Take the Fourier transform of each segment and
  track the Doppler shift by adding power in the
  frequency domain (Stack and Slide)

32
The Hough transform

• Robust pattern detection technique developed at
  CERN to look for patterns in bubble chamber
  pictures. Patented by IBM and used to detect
  patterns in digital images
• Look for patterns in the time-frequency plane
• Expected pattern depends on a,d, f0, fn

33
Hierarchical Hough transform strategy
Pre-processing
Divide the data set in N chunks
Template placing
Construct set of short FT (tSFT)
Candidates selection
Set upper-limit
Candidates selection
34
Incoherent Hough search Pipeline
Construct set of short FT (tSFTlt1800s)
Candidates selection
Set upper-limit
35
Peak selection in the t-f plane

• Input data Short Fourier Transforms (SFT) of
  time series
• For every SFT, select frequency bins i
  such exceeds some threshold rth
• ? time-frequency plane of zeros and ones
• p(rh, Sn) follows a ?2 distribution with 2
  degrees of freedom
• The false alarm and detection probabilities for a
  threshold rth are

36
Hough statistics

• After performing the HT using N SFTs, the
  probability that the pixel a,d, f0, fi has a
  number count n is given by a binomial
  distribution
• The Hough false alarm and false dismissal
  probabilities for a threshold nth
• ? Candidates selection
• For a given aH, the solution for nth is
• Optimal threshold for peak selection rth 1.6
  and a 0.20

37
The time-frequency pattern

• SFT data
• Demodulated data

Time at the SSB for a given sky position
38
Validation code-Signal only case- (f0500 Hz)

• bla

39
Validation code-Signal only case- (f0500 Hz)
40
Noise only case
41
Results on simulated data
Statistics of the Hough maps f0300Hz tobs60days
N1440 SFTs tSFT1800s a0.2231 ltngt
Na321.3074 ?n?(Na(1-a))15.7992
42
Number count probability distribution
1440 SFTs 1991 maps 58x58 pixels
43
Frequentist Analysis

• Perform the Hough transform for a set of points
  in parameter space la,d,f0,fi? S , given the
  data
• HT S ? N
• l ? n(l)
• Determine the maximum number count n
• n max (n(l)) l ? S
• Determine the probability distribution p(nh0)
  for a range of h0

44
Frequentist Analysis

• Perform the Hough transform for a set of points
  in parameter space la,d,f0,fi? S , given the
  data
• HT S ? N
• l ? n(l)
• Determine the maximum number count n
• n max (n(l)) l ? S
• Determine the probability distribution p(nh0)
  for a range of h0
• The 95 frequentist upper limit h095 is the
  value such that for repeated trials with a signal
  h0? h095, we would obtain n ? n more than 95
  of the time

• Compute p(nh0) via Monte Carlo signal
  injection, using l ? S , and ?0 ?0,2?, ?
  ?-?/4,?/4, cos ??-1,1.

45
Set of upper limit. Frequentist approach.
n395 ? ?0.295 ? h095
46
Comparison of sensitivity

• For a matched filter search directed at a single
  point in parameter space, smallest signal that
  can be detected with a false dismissal rate of
  10 and false alarm of 1 is
• For a Hough search with N segments, for same
  confidence levels and for large N
• For the Radon transform
• For, say N 2000, loss in sensitivity is only
  about a factor of 5 for a much smaller
  computational cost

47
Computational Engine

• Searchs offline at
• Medusa cluster (UWM)
• Merlin cluster (AEI)

48
Summary and future outlook

• S2 run (Feb 14, 2003 - Apr 14, 2003)
• Time-domain analysis of 28 known pulsars
• Broadband frequency-domain all-sky search
• ScoX-1 LMXB frequency-domain search
• Incoherent searches.
• S3 run (Oct 31, 2003 Jan 9, 2004)
• Time-domain analysis on more pulsars, including
  binaries
• Improved sensitivity LIGO/GEO run. Approaching
  spin-down limit for Crab pulsar
• Future
• Implement hierarchical analysis that layers
  coherent and incoherent methods
• Grid searches
• Einstein_at_home initiative for 2005 World Year of
  Physics