Transcript Searching for GW

Data analysis searches for
continuous gravitational waves
Alicia M. Sintes
Universitat de les Illes Balears
Pisa,1 June 2007
Content
• Gravitational waves from spinning neutron stars.
– Emission mechanisms
– Signal model
• Data analysis for continuous gravitational waves
– Bayesian analysis
– Frequentist approach
• Brief overview of the searches
– Directed pulsar search
– Wide-parameter search (all sky)
•
•
•
•
Coherent methods
Einstein@Home
Hierarchical strategies
Semi-coherent methods
• Summary of results and perspectives
2
VESF-school, Pisa, June 2007, A.M. Sintes
Rotating neutron stars
•
•
•
•
Neutron stars can form from the remnant of stellar collapse
Typical size of 10km, and are about 1.4 solar masses
Some of these stars are observed as pulsars
Gravitational waves from neutron stars could tell us about the equation of
state of dense nuclear matter
• Our galaxy might contain ~109 NS, of which ~105 are expected to be active
pulsars. Up to know ~1700 pulsars have been identified
• search for observed neutron stars
• all sky search (computing challenge)
(NASA/CXC/SAO)
3
NASA
VESF-school, Pisa, June 2007, A.M. Sintes
Neutron Stars Sources
• Great interest in detecting radiation:
physics of such stars is poorly
understood.
– After 40 years we still don’t
know what makes pulsars
pulse.
– Interior properties not
understood: equation of state,
superfluidity,
superconductivity, solid core,
source of magnetic field.
– May not even be neutron stars:
could be made of strange
matter!
4
VESF-school, Pisa, June 2007, A.M. Sintes
“Continuous” gravitational waves
from neutron stars
Various physical mechanisms could operate in a NS to produce interesting
levels of GW emission
As the signal-strength is generally expected to be weak, long integrations times
(several days to years) are required in order for the signal to be detectable in
the noise.
QuickTime™ and a
TIFF (Uncompressed) decompressor
are needed to see this picture.
Therefore this GW emission has to last for a long time, which
characterizes the class of ‘continuous-wave’ signals.
NS might also be interesting sources of burst-like GW emission,
f-mode or p-mode oscillations excited by a glitch
crustal torsion modes
5
VESF-school, Pisa, June 2007, A.M. Sintes
Emission mechanisms for continuous GW from
spinning NS in the LIGO-VIRGO frequency band
– Non-axisymmetric distortions
– Unstable oscillation modes in the
fluid part of the star
– Free precession
Low Mass X-Ray Binaries
Bumpy Neutron Star
Magnetic mountains
6
R-modes in accreting stars
Wobbling Neutron Star
VESF-school, Pisa, June 2007, A.M. Sintes
1) Non-axisymmetric distortions
A non-axisymmetric neutron star at distance a d, rotating with
frequency ν around the Izz axis emits monochromatic GWs of
frequency f = 2ν with an amplitude
2
4  2G Izz f gw
h0  4
c
d
Bumpy Neutron Star


Ixx  Iyy
Izz
•The strain amplitude h0 refers to a GW from an optimally
oriented source with respect to the detector

•The equatorial ellipticity is highly uncertain, ~10-7. In the
most speculative model can reach up to 10-4.
•Accreating neutron stars in binary systems can also have
large crust deformations
•Strong internal magnetic fields could produce
deformations of up to ~10-6. These deformations would
result in GW emission at f = ν and f = 2ν.
Magnetic mountains
7
VESF-school, Pisa, June 2007, A.M. Sintes
2) Non-axisymmetric instabilities
At birth or during accretion, rapidly rotating NS can be
subject to various non-axisymmetric instabilities, which
would lead to GW emission,
The r-mode instability has been proposed as a source of
GWs (with frequency f = 4ν/3) from newborn NS and
from rapidly accreting NS.
R-modes in accreting stars
3) Free precession
results in emission at (approximately) the rotation
rate ν and twice the rotation rate, i.e. f =ν+νprec
and f = 2ν.
QuickTime™ and a
TIFF (Uncompressed) decompressor
are needed to see this picture.
8
Wobbling
Neutron Star
VESF-school, Pisa, June 2007, A.M. Sintes
Loudest expected signal
From isolated neutron stars (Blandford’s argument):
Although there is great uncertainty in the physics of the GW emission mechanisms
and the strength of individual sources, one can argue for a statistical upper limit on
the expected strongest GW signals from the galactic population of neutron stars,
which is almost independent of individual source physics.
h0~410-24
Is the optimistic upper limit, that there is a 50% chance that the strongest signal
has at least that amplitude in the LIGO-VIRGO band, assuming NS are
uniformly distributed in the galactic disc and a constant overall galactic bith-rate.
Spindown limit for known pulsars
If the GW emission is powered by the rotational energy of the NS, then one obtains:
 ≤ sd, h0≤ hsd
QuickTime™ and a
TIFF (Uncompressed) decompressor
are needed to see this picture.
QuickTime™ and a
TIFF (Uncompressed) decompressor
are needed to see this picture.
9
VESF-school, Pisa, June 2007, A.M. Sintes
The signal from a NS
• The ‘periodic’ GW signal from a neutron
star:
(t)
h(t)  h0 A(t)e
• 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
T
•
(Signal-to-noise)2 ~

0
2
h (t)
dt
Sh (fgw )
10
VESF-school, Pisa, June 2007, A.M. Sintes
The expected signal at the detector
A gravitational wave signal we detect from a
NS will be:
– Frequency modulated by relative motion of
detector and source
– Amplitude modulated by the motion of the nonuniform antenna sensitivity pattern of the detector
R
11
VESF-school, Pisa, June 2007, A.M. Sintes
Signal received from
an isolated NS
h(t)  F (t;) h (t)  F (t; ) h  (t)
F+ and F are the strain antenna patterns. They depend on the orientation of the
detector and source and on the polarization of the waves.
 h   A cos(t)
h   A  sin(t)

f(n )
(t) 0  2 
(T(t)  T(t 0 ))n 1
n 0 (n  1)!
the phase of the received signal depends on the initial phase, the frequency
evolution 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, t the time at the
detector.
12
VESF-school, Pisa, June 2007, A.M. Sintes
Signal model:
isolated non-precessing NS
In the case of an isolated tri-axial neutron
star emitting at twice its rotational frequency
1
A   h0 (1 cos2  )
2
A   h0 cos
4  2G Izz f gw
h0  4
c
d
2
h0 - amplitude of the gravitational
wave signal
 - angle between the pulsar spin
axis and line of sight
Ixx  Iyy
- equatorial ellipticity

Izz
13
VESF-school, Pisa, June 2007, A.M. Sintes
Signals in noise.
The meaning of probability
The strain x(t) measured by a detector is mainly dominated by noise n(t), such
that even in the presence of a ignal h(t) we have
x(t)=n(t)+h(t)
Data analysis methods are not just simple recipes. We want tools capable of
—dealing with very faint sources
—handling very large data sets
—diagnosing systematic errors
—avoiding unnecessary assumptions
—estimating parameters and testing models
There are different ways of proceeding depending on the paradigm of statistics
used.
There are three popular interpretations of the word:
Probability as a measure of our degree of belief in a statement
Probability as a measure of limiting relative frequency of outcome of a set of identical
experiments
Probability as the fraction of favourable (equally likely) possibilities
We will call these the Bayesian, Frequentist and Combinatorial interpretations.
14
VESF-school, Pisa, June 2007, A.M. Sintes
Algebra of (Bayesian) probability
•
•
If there are two statements X and Y, then joint
probability
p(X,Y )  p(Y) p(X |Y)
where the vertical line denotes the conditional
statement “X given Y is true” – The Product Rule

…because
p(X,Y)=p(Y,X), we get
p(X |Y,I) 
p(X | I) p(Y | X,I)
p(Y | I)
Thomas Bayes
(1702 – 1761 AD)
which is called Bayes’ Theorem

p(model | data, I)
Prior
Likelihood
Posterior

p(data | model,I)  p(model | I)
p(data | I)
Bayes’ theorem is the
appropriate rule for updating
our degree of belief when we
have new data
Evidence

15
We can usually calculate all these terms
VESF-school, Pisa, June 2007, A.M. Sintes
Algebra of (Bayesian) probability
•
We can also deduce the marginal probabilities. If X and Y are propositions that can
take on values drawn from X  {x1, x 2,...,x n }, and Y  {y1, y 2,...,y m }, then
p(xi )  p(xi )

 p(y
j
| xi ) 
j1..m
 p(x )p(y
j1..m
i
j
| xi ) 
 p(x , y )
i
j
j1..m
=1

this gives use the probability of X when we don’t care about Y. In these
circumstances, Y is known as a nuisance parameter.
All these relationships can be smoothly extended from discrete probabilities to
probability densities, e.g.
•
p(x) 
 p(x,y)dy
where “p(y)dy” is the probability that y lies in the range y to y+dy.
16
VESF-school, Pisa, June 2007, A.M. Sintes
Bayesian parameter estimation
In the Bayesian approach, we can test our model
Posterior
Likelihood
Prior
p(model | data, I)  p(data | model,I)  p(model | I)
What we know now

Influence of our
observations
What we
knew before
• As our data improve (e.g. our sample increases), the posterior pdf narrows and
becomes less sensitive to our choice of prior.
• The posterior conveys our (evolving) degree of belief in different values of  ,
in the light of our data
• If we want to express our result as a single number we could perhaps adopt the
mean, median, or mode
• We can use the variance of the posterior pdf to assign an uncertainty for 
• It is very straightforward to define confidence intervals
17
VESF-school, Pisa, June 2007, A.M. Sintes
Bayesian parameter estimation
p( | data, I )
Bayesian confidence intervals
We are 95% sure that 
lies between 1 and 2
95% of area
under pdf
1
Note: the confidence interval is
not unique, but we can define
the shortest C.I.
2

18
VESF-school, Pisa, June 2007, A.M. Sintes
Frequentist framework
•
Recall that in frequentist (orthodox) statistics, probability is limiting relative
frequency of outcome, so :
only random variables can have frequentist probabilities
as only these show variation with repeated measurement. So we can’t talk about
the probability of a model parameter, or of a hypothesis. E.g., a measurement of a
mass is a random variable, but the mass itself is not.
•
So no orthodox probabilistic statement can be interpreted as directly referring to
the parameter in question! For example, orthodox confidence intervals do not
indicate the range in which we are confident the parameter value lies. That’s what
Bayesian intervals do.
19
VESF-school, Pisa, June 2007, A.M. Sintes
Frequentist hypothesis testing
– significance tests
The method goes like this:
– To test a hypothesis H1 consider another hypothesis, called the null hypothesis, H0, the
truth of which would deny H1. Then argue against H0…
– Use the data you have gathered to compute a test statistic obs which has a calculable
pdf if H0 is true. This can be calculated analytically or by Monte Carlo methods.
– Look where your observed value of the statistic lies in the pdf, and reject H0 based on
how far in the wings of the distribution you have fallen.
p(|H0)
Reject H0 if your result lies in here

20
VESF-school, Pisa, June 2007, A.M. Sintes
Frequentist hypothesis testing
– significance tests
•
H0 is rejected at the x% level if x% of the probability lies to the right of the
observed value of the statistic (or is ‘worse’ in some other sense):
p(|H0)
X% of the
area

obs
•
and makes no reference to how improbable the value is under any alternative
hypothesis (not even H1!).
We will chose a critical region of *, so that if obs>* we reject H0 and
therefore accept H1.
p(|H0)
21
p(|H1)


VESF-school, Pisa, June 2007, A.M. Sintes
Frequentist hypothesis testing
QuickTime™ and a
TIFF (Uncompressed) decompressor
are needed to see this picture.
QuickTime™ and a
TIFF (Uncompressed) decompressor
are needed to see this picture.
Type II error
Type I error

–
–
A Type I error occurs when we reject the null hypothesis when it is true (false alarm)
A Type II error occurs when we accept the null hypothesis when it is false (false dismissal)
both of which we should strive to minimise. According to the Neyman-Pearson lemma, this optimal test is
the so-called likelihood ratio:
QuickTime™ and a
TIFF (Uncompressed) decompressor
are needed to see this picture.
•
For Gaussian noise, one finds
QuickTime™ and a
TIFF (Uncompressed) decompressor
are needed to see this picture.
QuickTime™ and a
TIFF (Uncompressed) decompressor
are needed to see this picture.
which is the well-known expression for the matched-filtering amplitude. If some of the parameters of the signal
h(t;A,) are unknown, one has to find the maximum of lnΛ as a function of the unknown parameters
22
VESF-school, Pisa, June 2007, A.M. Sintes
Calibrated output: LIGO noise history
Integration times
S1 - L1 5.7 days, H1 8.7
days, H2 8.9 days
S2 - L1 14.3 days, H1 37.9
days, H2 28.8 days
S3 - L1 13.4 days, H1 45.5
days, H2 42.1 days
S4 - L1 17.1 days, H1 19.4
days, H2 22.5 days
S5 (so far...) >1 year
23
VESF-school, Pisa, June 2007, A.M. Sintes
Calibrated output: GEO noise history
24
VESF-school, Pisa, June 2007, A.M. Sintes
Continuous wave searches
•
Signal parameters: position (may be known), inclination angle, [orbital parameters in case
of a NS in a binary system], polarization, amplitude, frequency (may be known), frequency
derivative(s) (may be known), initial phase.
•
Most sensitive method: coherently correlate the data with the expected signal (template)
and inverse weights with the noise. If the signal were monochromatic this would be
equivalent to a FT.
– Templates: we assume various sets of unknown parameters and correlate the data against these
different wave-forms.
– Good news: we do not have to search explicitly over polarization, inclination, initial phase and
amplitude.
•
Because of the antenna pattern, we are sensitive to all the sky. Our data stream has signals
from all over the sky all at once. However: low signal-to-noise is expected. Hence
confusion from many sources overlapping on each other is not a concern.
•
Input data to our analyses:
– A calibrated data stream which with a better than 10% accuracy, is a measure of the GW
excitation of the detector. Sampling rate 16kHz (LIGO-GEO, 20kHz VIRGO), but since the high
sensitivity range is 40-1500 Hz we can downsample at 3 kHz.
25
VESF-school, Pisa, June 2007, A.M. Sintes
Four neutron star populations
and searches
•
Known pulsars
• Position & frequency evolution known (including derivatives, timing noise, glitches, orbit)
•
Unknown neutron stars
• Nothing known, search over position, frequency & its derivatives
•
Accreting neutron stars in low-mass x-ray binaries
• Position known, sometimes orbit & frequency
•
Known, isolated, non-pulsing neutron stars
• Position known, search over frequency & derivatives
• What searches?
– Targeted searches for signals from known pulsars
– Blind searches of previously unknown objects
– Coherent methods (require accurate prediction of the phase evolution of the signal)
– Semi-coherent methods (require prediction of the frequency evolution of the signal)
What drives the choice? The computational expense of the search
26
VESF-school, Pisa, June 2007, A.M. Sintes
Coherent detection methods
There are essentially two types of coherent searches that are performed
•
Frequency domain
Time domain
Conceived as a module in a hierarchical search
process signal to remove frequency variations due to
Earth’s motion around Sun and spindown
Matched filtering techniques.
Aimed at computing a detection
statistic.
These methods have been implemented in
the frequency domain (although this is not
necessary) and are very computationally
efficient.
•
Best suited for large parameter
space searches
•
Standard Bayesian analysis,
as fast numerically but provides natural
parameter estimation
•
Best suited to target known
objects, even if phase evolution
is complicated
•
Efficiently handless missing data
•
Upper limits interpretation:
Bayesian approach
(when signal characteristics are uncertain)
•
Frequentist approach used to
cast upper limits.
27
VESF-school, Pisa, June 2007, A.M. Sintes
Summary of directed
pulsar searches
•
•
Within the LIGO sensitive band (gw > 50 Hz) there are currently 163
known pulsars
We have rotational and positional parameter information for 97 of these
• S1 (LIGO and GEO: separate analyses)
– Upper limit set for GWs from J1939+2134 (h0<1.4 x 10-22)
– Phys. Rev. D 69, 082004 (2004)
• S2 science run (LIGO: 3 interferometers coherently, TDS)
– End-to-end validation with 2 hardware injections
– Upper limits set for GWs from 28 known isolated pulsars
– Phys. Rev. Lett. 94, 181103 (2005)
• S3 & S4 science runs (LIGO and GEO: up to 4 interferometers
coherently, TDS)
– Additional hardware injections in both GEO and LIGO
– Add known binary pulsars to targeted search
– Full results with total of 93 (33 isolated, 60 binary) pulsars
• S5 science run (ongoing, TDS)
– 32 known isolated, 41 in binaries, 29 in globular clusters
28
VESF-school, Pisa, June 2007, A.M. Sintes
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
(h|h). Fixing ,  and 0 , for every h0, we can obtain a pdf curve:
p(2F|h0)
• 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*
C(h0 ) 


2F*
p(2F | h0 )d(2F)
• Use signal injection Monte Carlos to measure Probability Distribution
Function (PDF) of F

29
VESF-school, Pisa, June 2007, A.M. Sintes
Measured PDFs for the F statistic
with fake injected worst-case signals at nearby frequencies
h0 = 1.9E-21
h0 = 2.7E-22
S1
Note:
hundreds
of
thousands
of injections
were
needed to
get such
nice clean
statistics!
95%
95%
2F* = 1.5: Chance probability 83%
2F*
2F*
2F* = 3.6: Chance probability 46%
h0 = 5.4E-22
95%
h0 = 4.0E-22
95%
2F*
2F* = 6.0: chance probability 20%
2F*
2F* = 3.4: chance probability 49%
30
VESF-school, Pisa, June 2007, A.M. Sintes
F-statistics
QuickTime™ and a
TIFF (Uncompressed) decompressor
are needed to see this picture.
QuickTime™ and a
TIFF (Uncompressed) decompressor
are needed to see this picture.
We can express h(t) in terms
of amplitude A {A+, A,
, 0} and Doppler
parameters 
QuickTime™ and a
TIFF (Uncompressed) decompressor
are needed to see this picture.
QuickTime™ and a
TIFF (Uncompressed) decompressor
are needed to see this picture.
QuickTime™ and a
TIFF (Uncompressed) decompressor
are needed to see this picture.
31
VESF-school, Pisa, June 2007, A.M. Sintes
F-Statistics
• Analytically maximize the
likelihood over A
QuickTime™ and a
TIFF (Uncompressed) decompressor
are needed to see this picture.
QuickTime™ and a
TIFF (Uncompressed) decompressor
are needed to see this picture.
QuickTime™ and a
TIFF (Uncompressed) decompressor
are needed to see this picture.
QuickTime™ and a
TIFF (Uncompressed) decompressor
are needed to see this picture.
QuickTi me™ and a
TIFF (Uncompressed) decompressor
are needed to see thi s pi cture.
QuickTime™ and a
TIFF (Uncompressed) decompressor
are needed to see this picture.
32
VESF-school, Pisa, June 2007, A.M. Sintes
Time domain target search
• 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
y(t;a)  14 h0F (t,)(1 cos2  )e2i 0  12 ih0F (t,)cose 2i 0

33
VESF-school, Pisa, June 2007, A.M. Sintes
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.
p(a |{Bk })  p(a)  p({Bk } | a)
posterior
• Thelikelihood 
prior
exp(-c2 /2),
likelihood
where
c 2 (a)  
k
Bk  y(t;a)
2
k
• To get the posterior PDF for h0, marginalizing with respect to the
nuisance parameters cos ,  and 0 given the data Bk
p(h0 |{Bk })   e

c 2 / 2
d0ddcos
34
VESF-school, Pisa, June 2007, A.M. Sintes
Posterior PDFs for
CW time domain analyses
Simulated injection
at 2.2 x10-21
p
shaded area =
95% of total area
p
35
VESF-school, Pisa, June 2007, A.M. Sintes
Progression of targeted pulsars
upper limits
•
•
Amplitudes of < 10-25 and
ellipticities <10-6 for many
objects
Our most stringent ellipticities
(10-7) are starting to reach into
the range of neutron star
structures for some neutronproton-electron models
h0
r 1Hz2 1038 kgm2
  0.237 24
10 1kpc  2
Izz
Crab pulsar
36
VESF-school, Pisa, June 2007, A.M. Sintes
Preliminary S5 results
Used parameters provided by Pulsar Group, Jodrell Bank
Joint 95% upper limits from first ~13 months
of S5 using H1, H2 and L1 (97 pulsars)
Lowest h0 upper limit:
PSR J1435-6100 (ngw = 214.0 Hz, r = 3.3 kpc) h0_min = 4.2x10-26
Lowest ellipticity upper limit:
PSR J2124-3358 (ngw = 405.6Hz, r = 0.25 kpc)  = 9.6x10-8
37
APS Meeting, Jacksonville, April
2007
Due to pulsar glitches the Crab
pulsar result uses data up to
the glitch on 23 Aug 2006,
and the PSRJ0537-6910 result
uses only three months of data
between two glitches on 5th
May and 4th Aug 2006
VESF-school, Pisa, June 2007, A.M. Sintes
preliminary
38
APS Meeting, Jacksonville, April
2007
• Black curve
represents one
full year of data
for all three
interferometers
running at design
sensitivity
• Blue stars
represent pulsars
for which we are
reasonably
confident of
having phase
coherence with
the signal model
• Green stars
represent pulsars
for which there is
uncertainty about
phase coherence
VESF-school, Pisa, June 2007, A.M. Sintes
Blind searches and coherent
detection methods
• Coherent methods are the most sensitive methods (amplitude
SNR increases with sqrt of observation time) but they are the
most computationally expensive,
why?
– Our templates are constructed based on different values of the signal
parameters (e.g. position, frequency and spindown)
– The parameter resolution increases with longer observations
– Sensitivity also increases with longer observations
– As one increases the sensitivity of the search, one also increases the
number of templates one needs to use.
39
VESF-school, Pisa, June 2007, A.M. Sintes
Number of templates
The number of templates grows dramatically with the
coherent integration time baseline and the computational
requirements become prohibitive
[Brady et al., Phys.Rev.D57 (1998)2101]
40
VESF-school, Pisa, June 2007, A.M. Sintes
S2 run: Coherent search for
unknown isolated sources and Sco-X1
•
•
•
•
•
Entire sky search
Fully coherent matched filtering
160 to 728.8 Hz
df/dt < 4 x 10-10 Hz/s
10 hours of S2 data;
computationally intensive
• 95% confidence upper limit on the
GW strain amplitude range from
6.6x10-23 to 1.0x10-21 across the
frequency band
•
•
•
•
•
•
Scorpius X-1
Fully coherent matched filtering
464 to 484 Hz, 604 to 624 Hz
df/dt < 1 x 10-9 Hz/s
6 hours of S2 data
95% confidence upper limit on the
GW strain amplitude range from
1.7x10-22 to 1.3x10-21 across the
two 20 Hz wide frequency bands
• See gr-qc/0605028
41
VESF-school, Pisa, June 2007, A.M. Sintes
Computational Engine
Searchs offline at:
• Medusa, Nemo clusters
(UWM)
• Merlin, Morgane cluster
(AEI)
• Tsunami (B’ham)
• Others
42
VESF-school, Pisa, June 2007, A.M. Sintes
Einstein@home
•
•
•
•
•
•
Like SETI@home, but for
LIGO/GEO data
American Physical Society
(APS) publicized as part of
World Year of Physics (WYP)
2005 activities
Use infrastructure/help from
SETI@home developers for
the distributed computing
parts (BOINC)
Goal: pulsar searches using
~1 million clients. Support
for Windows, Mac OSX,
Linux clients
From our own clusters we can
get ~ thousands of CPUs.
From Einstein@home hope to
get order(s) of magnitude
more at low cost
Currently : ~140,000 active
users corresponding to about
80Tflops
http://einstein.phys.uwm.edu/
43
VESF-school, Pisa, June 2007, A.M. Sintes
Coherent searches on
Einstein@home
• Public distributed computing project to look for isolated
pulsars in LIGO/GEO data ~ 80 TFlops 24/7
• Initially making use of coherent F-statistic method
S3 - no spindown
• Using segment lengths of 10 hours
• No evidence of strong pulsar signals
• Outliers are consistent with instrumental artifacts or bad
bands. None of the low significance remaining candidates
showed up in follow-up on S4 data.
S4 - one spindown parameter, up to f/fdot ~ 10,000 yr
• Using segment lengths of 30 hours
• Analysis took ~ 6 months
• Post-processing complete. Currently under review.
S5 (R1)- similar to S4
• Faster more efficient application
• Currently in post-processing stage
44
VESF-school, Pisa, June 2007, A.M. Sintes
User/Credit History
45
http://www.boincsynergy.com/stats/
VESF-school, Pisa, June 2007, A.M. Sintes
http://www.boincstats.com/
Current performance
46
Einstein@Home is currently getting 84 Tflops
VESF-school, Pisa, June 2007, A.M. Sintes
All-Sky surveys for unknown
gravity-wave emitting pulsars
It is necessary to search for every signal template distinguishable in parameter
space. Number of parameter points required for a coherent T=107s search
[Brady et al., Phys.Rev.D57 (1998)2101]:
Class
f (Hz)
t(Yrs)
Ns
Directed
All-sky
Slow-old
<200
>103
1
3.7x106
1.1x1010
Fast-old
<103
>103
1
1.2x108
1.3x1016
Slow-young
<200
>40
3
8.5x1012
1.7x1018
Fast-young
<103
>40
3
1.4x1015
8x1021
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.
47
VESF-school, Pisa, June 2007, A.M. Sintes
Coherent wide-parameter searches
• The second effect of the large number of templates Np is to reduce the
sensitivity compared to a targeted search with the same observation time
and false-alarm probability: increasing the number of templates increases
the number of expected false-alarm candidates at fixed detection
threshold. Therefore the detection-threshold needs to be raised to maintain
the same false-alarm rate, thereby decreasing the sensitivity.
• Note that increasing the number of equal-sensitivity detectors N improves
the SNR in the same way as increasing the integration time Tobs. However,
increasing the number of detectors N does — contrary to the observation
time Tobs — not increase the required number of templates Np, which
makes this the computationally cheapest way to improve the SNR of
coherent wide-parameter searches.
48
VESF-school, Pisa, June 2007, A.M. Sintes
Hierarchical strategies
Pre-processing
raw data
GEO/LIGO
Divide the data set in N
chunks
Template
placing
Coherent search (a,d,fi)
Construct set of short
FT (tSFT)
in a frequency band
Incoherent search
Peak selection
in t-f plane
Set upper-limit
Hough transform
(a,d,f0, fi)
Candidates
selection
Candidates
selection
49
VESF-school, Pisa, June 2007, A.M. Sintes
Hierarchical method for ‘blind’ searches
h-reconstructed data
Data quality
Data quality
SFDB
SFDB
Average spectrum
estimation
Average spectrum
estimation
peak map
peak map
Hough transf.
Hough transf.
candidates
coincidences
The procedure involves two or
more data sets belonging to a
single or more detectors
candidates
coherent
step
50
events
VESF-school, Pisa, June 2007, A.M. Sintes
Incoherent power-sum methods
•
The idea is to perform a search over the total observation
time using an incoherent (sub-optimal) method:
•
Three methods have been developed to search for
cumulative excess power from a hypothetical periodic
gravitational wave signal by examining successive spectral
estimates:
– Stack-slide (Radon transform)
– Hough transform
– Power-flux method
They are all based on breaking up the data into segments,
FFT each, producing Short (30 min) Fourier Transforms
(SFTs) from h(t), as a coherent step (although other
coherent integrations can be used if one increasing the
length of the segments), and then track the frequency drifts
due to Doppler modulations and df/dt as the incoherent
step.
Frequency
Tcoh
Time
51
VESF-school, Pisa, June 2007, A.M. Sintes
Differences among the
incoherent methods
What is exactly summed?
• StackSlide – Normalized power (power divided by
estimated noise)
 Averaging gives expectation of 1.0 in absence of signal
• Hough – Weighted binary counts (0/1 = normalized power
below/above SNR), with weighting based on antenna
pattern and detector noise
• PowerFlux – Average strain power with weighting based on
antenna pattern and detector noise
 Signal estimator is direct excess strain noise
(circular polarization and 4 linear polarization projections)
52
VESF-school, Pisa, June 2007, A.M. Sintes
Peak selection in the t-f plane
• Input data: Short Fourier Transforms (SFT) of time series
•
For every SFT, select frequency bins i such x˜(f)  n˜ (f)  h˜(f)
| x˜(fi ) |2
ri 
| n˜ (fi ) |2
exceeds some threshold rth
time-frequency plane of zeros and ones


• p(r|h, Sn) follows a c2 distribution with 2 degrees of freedom:
4 | h˜(fi ) |2
i 
Sn (fi )TSFT
i
ri  1 
2
 r 2  1  i
• The false alarm and detection probabilities for a threshold rth are

a ( r th | Sn )
 p(r | 0,S
n
)dr  e r th ,
r th
( r th | h,Sn )
53


 p(r | h,S
n
)dr
r th
VESF-school, Pisa, June 2007, A.M. Sintes
Standard 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:
N  n
P(n | p,N)  p (1 p) N  n
n 
n  Np,  n  Np(1  p)
2
a signal absent
p  
 signal present
• The Hough false alarm and false dismissal probabilities for a threshold nth

N  n
n
a H (n th ,a,N)   a (1 a ) N
n
n n th  
N
 H (n th ,,N)
N  n
n  (1 )N n
n 0
n th 1
 Candidates selection
• For a given aH, the solution for nth is
nth  Na  2Na(1 a) erfc1(2a H )

• Optimal threshold for peak selection : rth ≈1.6 and a≈0.20
54
VESF-school, Pisa, June 2007, A.M. Sintes
Noise only case
55
VESF-school, Pisa, June 2007, A.M. Sintes
Frequentist upper limit
•
Perform the Hough transform for a set of points in
parameter space ={a,d,f0,fi} S , given the data:
HT: S  N
  n()
Determine the maximum number count n*
n* = max (n()):  S
Determine the probability distribution p(n|h0) 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
N
0.95
 pn|h )
95%
0
n n*
Compute p(n|h0) via Monte Carlo signal injection, using
  S , and 0 [0,2],  [-/4,/4], cos [-1,1].
56
VESF-school, Pisa, June 2007, A.M. Sintes
Number count distribution
for signal injections
L1: 200-201 Hz, n* =202
1000 injections
p(n|h0) ideally binomial
for a target search, but:
0.1%
30.5%
87.0%
1
– Amplitude modulation of
the signal for different SFTs
– Different sensitivity for
different sky locations
– Random mismatch between
signal & templates
‘smear’ out the binomial
distributions
57
VESF-school, Pisa, June 2007, A.M. Sintes
Hough S2: UL Summary
Feb.14-Apr.14,2003
•
•
S2 analysis covered 200-400Hz, over the
whole sky, and 11 values of the first
spindown (Δf = 5.55×10 – 4 Hz, Δf1 = –
1.1×10– 10 Hz s– 1)
Templates: Number of sky point templates
scales like (frequency)2
–
–
–
•
•
•
1.5×105 sky locations @ 300 Hz
1.9×109 @ 200-201 Hz
7.5×109 @ 399-400 Hz
Three IFOs analyzed separately
No signal detected
Upper limits obtained for each 1 Hz band
by signal injections: Population-based
frequentist limits on h0 averaging over sky
location and pulsar orientation
58
Detector
Frequency (Hz)
h095%
L1
200-201
4.43x10-23
H1
259-260
4.88x10-23
H2
258-259
8.32x10-23
VESF-school, Pisa, June 2007, A.M. Sintes
Improvements for S4
• We use a weighted Hough to give more weight to SFTs having greater SNR.
Weights are proportional to the beam pattern functions and inversely
proportional to the SFT noise floor. 2
2
F )  F )

w 
(i)

i
(i)

Sn(i)
• Number count n is not an integer anymore
N
n   wi ni
i1
N
w
i
N
i1
• Using the weights does not lead to any loss in computational efficiency or
 also been generalized to the Multi-IFO case
robustness. It has
• Mean number count is unchanged due to normalization of weights. Standard
deviation always increases:

n  Na  N exp rth )

 n  w a (1  a )
• Number count threshold for a given false alarm:
2
nth  Na  2 w a (1  a ) erfc1 (2a H )
59
VESF-school, Pisa, June 2007, A.M. Sintes
Contribution of different IFOs
Figure plots the relative noise weights from H1, H2 and L1


IFO
wi
All
wi
60
VESF-school, Pisa, June 2007, A.M. Sintes
Loudest events for every 0.25 Hz
Multi interferometer case 50-1000 Hz
Significance defined as s=(nmax-<n>)/
61
VESF-school, Pisa, June 2007, A.M. Sintes
The S4 Hough search
•
•
•
•
•
•
•
As before, input data is a set of N
1800s SFTs (no demodulations)
Weights allow us to use SFTs from
all three IFOs together:
1004 SFTS from H1, 1063 from
H2 and 899 from L1
Search frequency band 50-1000Hz
1 spin-down parameter. Spindown
range [-2.2,0]×10-9 Hz/s with a
resolution of 2.2×10-10 Hz/s
All sky search
All-sky upper limits set in 0.25 Hz
bands
Multi-IFO and single IFOs have
been analyzed
Best UL
for L1: 5.9×10-24
for H1: 5.0×10-24
for Multi H1-H2-L1: 4.3×10-24
62
VESF-school, Pisa, June 2007, A.M. Sintes
Improvements due to the weights
Comparison of the All-sky 95%
upper limits obtained by MonteCarlo injections for the multiIFO case.
The average improvement by
using weights in this band is
9.25% for the multi-IFO case,
but only ~6% for the single IFO
63
VESF-school, Pisa, June 2007, A.M. Sintes
StackSlide/Hough/PowerFlux
differences
StackSlide
Hough
PowerFlux
Windowing
Tukey
Tukey
Hann
Noise estimation
Median-based floor Median-based floor Time/frequency
tracking
tracking
decomposition
Line handling
Cleaning
Cleaning
Skyband exclusion
Antenna pattern
weighting
No
Yes
Yes
Noise weighting
No
Yes
Yes
Spindown step size
2 x 10-10 Hz/s
2 x 10-10 Hz/s
Freq dependent
Limit at every
skypoint
No
No
Yes
Upper limit type
Population-based
Population-based
Strict frequentist
64
VESF-school, Pisa, June 2007, A.M. Sintes
H1 (Hanford 4-km) and Multi-IFO Results
PRELIMINARY
65
VESF-school, Pisa, June 2007, A.M. Sintes
L1 (Livingston 4-km) and Multi-IFO Results
PRELIMINARY
66
VESF-school, Pisa, June 2007, A.M. Sintes
S5 incoherent searches
preliminary PowerFlux results
67
VESF-school, Pisa, June 2007, A.M. Sintes
E@H-Hierarchical search
Basic algorithm
QuickTime™ and a
TIFF (Uncompressed) decompressor
are needed to see this picture.
• Started test-run (~3 months)
– currently working on improving stability & reliability
• Will be followed by a full 1-year run
• Should permit a search that extends hundreds of pc into the Galaxy
• This should become the most sensitive blind CW search possible with
current knowledge and technology
68
VESF-school, Pisa, June 2007, A.M. Sintes
Sensitivity estimate of different
E@H Runs
QuickTime™ and a
TIFF (Uncompressed) decompressor
are needed to see this picture.
69
VESF-school, Pisa, June 2007, A.M. Sintes
Searches for Continuous Waves,
present, past and future
QuickTime™ and a
TIFF (Uncompressed) decompressor
are needed to see this picture.
70
VESF-school, Pisa, June 2007, A.M. Sintes
Conclusions
• Analysis of LIGO data is in full swing, and results from
LIGO searches from science runs 4, 5 are now appearing.
– Significant improvements in interferometer sensitivity since S3.
– In the process of accumulating 1 year of data (S5).
– Known pulsar searches are beginning to place interesting upper limits
in S5
– All sky searches are under way and exploring large area of parameter
space
• Veto strategies and data conditioning are becoming more relevant in
order to reduce the selection threshold and improve sensitivity
71
VESF-school, Pisa, June 2007, A.M. Sintes