F-statistic Search for white-dwarf binaries - Caltech (Tapir)

Mock LISA Data Challenge 2: F-statistic Search for white-dwarf binaries Note for results submission

R. Prix and J. T. Whelan June 16, 2007

Contents 1 General search method 1.1 Multi-detector F-statistic . . . . . . . . . . . . . . . . . . . . 1.2 Parameter-estimation . . . . . . . . . . . . . . . . . . . . . . .

1 2 4

2 Application to LISA and our MLDC pipeline 2.1 MLDC conventions for amplitude parameters . 2.2 TDI and long-wavelength approximation . . . 2.3 Wide-parameter search grid . . . . . . . . . . 2.4 MLDC2 WDB Pipeline . . . . . . . . . . . . .

6 6 6 8 9

1

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

General search method

We use the F-statistic, a coherent matched-filtering detection statistic first introduced by Jaranowski et al. [6] in the context of the search for continuouswave signals in ground based detectors. This method has been implemented by the LIGO Scientific collaboration in LAL/LALApps [8], and is currently used in the search for quasi-periodic GW signals from spinning neutron stars (e.g. see [1]). The generalization of the F-statistic to a coherent multidetector search was first obtained by Cutler and Schutz [4]. The application of the F-statistic to the search of continuous-wave sources (such as galactic white-dwarf binaries) using LISA was first discussed in Królak et al. [7]. The multi-detector F-statistic has been implemented in LALapps, in the code ComputeFStatistic v2, which we are using for the present analysis.

1

This is the same method and code that has been used previously by us on the first round of the MLDC [3, 9].

1.1

Multi-detector F-statistic

As shown in [6], the dimensionless strain signal sX (t) of a continuous gravitational wave at detector X can be represented in the form sX (t) =

4 X

Aµ hX µ (t) ,

(1.1.1)

µ=1

in terms of four signal-amplitudes Aµ , which are independent of the detector X, and the detector-dependent basis waveforms hX µ (t). The four amplitudes µ A can be expressed in terms of two polarization amplitudes A+ , A× , the initial phase φ0 in the solar-system barycenter (SSB) at a reference time τref , and the polarization angle ψ of the wave frame with respect to the equatorial coordinate system, namely A1 A2 A3 A4

= A+ cos φ0 cos 2ψ − A× sin φ0 sin 2ψ , = A+ cos φ0 sin 2ψ + A× sin φ0 cos 2ψ , = −A+ sin φ0 cos 2ψ − A× cos φ0 sin 2ψ , = −A+ sin φ0 sin 2ψ + A× cos φ0 cos 2ψ .

(1.1.2a) (1.1.2b) (1.1.2c) (1.1.2d)

We can further relate the two polarization amplitudes A+ and A× to the overall amplitude h0 and the inclination angle ι of the quadrupole rotation axis with respect to the line of sight, namely 1 A+ = h0 1 + cos2 ι , 2

A× = h0 cos ι .

(1.1.3)

The four basis waveforms hX µ (t) can be written as X X X X X hX 1 (t) = a (t) cos φ (t) , h2 (t) = b (t) cos φ (t) , X X X X X hX 3 (t) = a (t) sin φ (t) , h4 (t) = b (t) sin φ (t) ,

(1.1.4)

where aX (t) and bX (t) are the antenna-pattern functions (see Eqs.(12,13) of [6]), and φX (t) is the signal phase at the detector X. The antenna-pattern functions aX (t), bX (t) depend on the sky position of the GW source (which is equivalent to the propagation direction b k of the wave), and on the location and orientation of the detector X. The phase φX (t) also depends on the intrinsic phase parameters, ω say, of the signal. In the case of continuous 2

waves from isolated neutron stars, ω would only consist of the s + 1 spin parameters, i.e. ω = {f (k) }sk=0 , where f (k) is the k-th time-derivative of the intrinsic signal frequency in the SSB. In the following we denote the set of “Doppler parameters” (i.e. the parameters affecting the time evolution of the phase) by λ ≡ {b k, ω}, as opposed to µ µ the four “amplitude parameters” {A} = A . Using the multi-detector notation of [4, 7], we write vectors in “detectorspace” in boldface, i.e. {s}X = sX , and so the signal model (1.1.1) can be written as s(t; A, λ) = Aµ hµ (t; λ) , (1.1.5) with implicit summation over repeated amplitude indices, µ ∈ {1, 2, 3, 4}. The multi-detector scalar product is defined as Z ∞ −1 (x|y) ≡ x eX (f ) SXY (f ) yeY∗ (f ) df , (1.1.6) −∞

where x e(f ) denotes the Fourier transform of x(t). We use implicit summation over repeated detector indices, and the inverse noise matrix is defined by −1 YZ S = δXZ . In the case of uncorrelated noise, where S XY = S X δ XY , the SXY scalar product simplifies to X (xX |y X ) , (1.1.7) (x|y) = X

in terms of the usual single-detector scalar product Z ∞ X x e (f ) yeX∗ (f ) X X (x |y ) ≡ df . S X (f ) −∞

(1.1.8)

With the signal model (1.1.1), the log-likelihood ratio is found as 1 ln Λ(x; A, λ) = Aµ xµ − Aµ Mµν Aν , 2

(1.1.9)

xµ (λ) ≡ (x|hµ ) , Mµν (λ) ≡ (hµ |hν ) .

(1.1.10) (1.1.11)

where we defined

We see that the likelihood ratio (1.1.9) can be maximized analytically with respect to the unknown amplitudes Aµ , resulting in the maximum-likelihood estimators (1.1.12) AµMLE = Mµν xν . 3

Substituting this into the detection statistic, we obtain the so-called Fstatistic, namely 2F(x; λ) ≡ xµ Mµν xν , (1.1.13) where Mµν ≡ {M−1 }µν , i.e. Mµα Mαν = δµν . Let us consider the case where the target Doppler parameters λ are perfectly matched to the signal λs , we find the expectation value of the F-statistic as E[2F] = 4 + SNR2 ,

(1.1.14)

in terms of the “optimal” signal-to-noise ratio SNR, which is expressible as SNR2 = sµ Mµν sν = Aµ Mµν Aν = (s|s) .

1.2

(1.1.15)

Parameter-estimation

From the expression (1.1.12) for the maximum-likelihood amplitudes Aµ in terms of the measured Fa , Fb , we can infer the signal-parameters A+ , A× (or equivalently h0 , cos ι) and ψ, φ0 , by using (1.1.3) and (1.1.2). We compute the two quantities A2s ≡

4 X

(Aµ )2 = A2+ + A2× ,

µ=1 1 4

Da ≡ A A − A2 A3 = A+ A× , which can easily be solved for A+ , A× , namely p 2A2+,× = A2s ± A4s − 4Da2 ,

(1.2.1) (1.2.2)

(1.2.3)

where our convention here is |A+ | ≥ |A× |, cf. (1.1.3), and therefore the ’+’ solution is A+ , and the 0 −0 is A× . The sign of A+ is always positive by this convention, while the sign of A× is given by the sign of Da , as can be seen from (1.2.2). Note that the discriminant in (1.2.3) is also expressible as p disc ≡ A4s − 4 Da2 = A2+ − A2× ≥ 0 . (1.2.4) Having computed A+ , A× , we can now also obtain ψ and φ0 , namely defining β ≡ A× /A+ , and b1 ≡ A4 − βA1 , b2 ≡ A3 + βA2 , b3 ≡ βA4 − A1 , 4

(1.2.5) (1.2.6) (1.2.7)

we easily find 1 ψ = atan 2

b1 b2

.

(1.2.8)

.

(1.2.9)

and φ0 = atan

b2 b3

The amplitudes Aµ are seen from (1.1.2) to be invariant under the following gauge-transformation, namely simultaneously {ψ → π/2, φ0 → φ0 + π}. Applying this twice, and taking account of the trivial gauge-freedom by 2π, this also contains the invariance ψ → ψ +π. Note that there is still an overall sign-ambiguity in the amplitudes Aµ , which can be determined as follows: compute a ’reconstructed’ A1r from (1.1.2) using the estimates A+,× and ψ, φ0 , and compare its sign to the original estimate A1 of (1.1.12). If the sign differs, the correct solution is simply found by replacing φ0 → φ0 + π. In order to fix a unique gauge, we restrict the quadrant of ψ to be ψ ∈ [−π/4, π/4) (in accord with the TDS convention), which can always be achieved by the above gauge-transformation, while φ0 remains unconstrained in φ0 ∈ [0, 2π). Converting A+ , A× into h0 and µ ≡ cos ι is done by solving (1.1.3), which yields q A2+ − A2× ,

h0 = A+ +

(1.2.10)

where we only kept the ’+’ solution, as we must have h0 > A+ . Finally, µ = cos ι is simply given by cos ι = A× /h0 . We know that the errors dxµ satisfy (assuming Gaussian noise): E[dxµ dxν ] = Mµν .

(1.2.11)

As a consequence of (1.1.12), we therefore obtain the covariance-matrix of the estimation-errors dAµ as E[dAµ dAν ] = Mµν ,

(1.2.12)

which corresponds to the Cramér-Rao bound, and Mµν is seen to be the inverse Fisher-matrix. The corresponding Fisher matrix for the variables {h0 , cos ι, ψ, φ0 } is simply obtained from the above together with the appropriate Jacobian accounting for the variables-transformation from Aµ .

5

2

Application to LISA and our MLDC pipeline

2.1

MLDC conventions for amplitude parameters

Unfortunately, the MLDC conventions for the amplitude parameters differ from the above standard LIGO/CW definitions for {h0 , cos ι, ψ, φ0 }. Here we only summarize without derivation how the “translation” is performed: • MLDC “Amplitude” = h0 /2 • MLDC “Inclination” = π − ι • MLDC “Polarization” = π/2 − ψ • MLDC “InitialPhase” = φ0 Note: There still seems to be an overall sign-difference (corresponding to a difference of π in the “InitialPhase”) between our estimated amplitudeparameters and those quoted in the MLDC keys, which is not understood. This sign-differences was already noted in MLDC1 and is still present in MLDC2, and we therefore simply “correct” our quoted “InitialPhase” by π.

2.2

TDI and long-wavelength approximation 2 v

-

...... ....... .. ....... .... ...... ... ....... ... ....... ... ....... . . . . . . ... .... ....... ... ...... ... ....... . . . . . . ... .. ....... . ... . . . . . ... ... ....... . . . . ... . .. ....... .. . . . . . . ... .. ....... . ... . . . . . ... ... ....... . . . . ... . .. ... 3 ..................... ... . . . . . ... . . ... ... ....... . . . . . .. ... 1 ....... . . ... . . . . .. ... ....... . . . . . . ... ... ....... .. . . . . . ... .. ....... . . . . . . . . ... ...... ... ....... ....... ... ....... ... ...... ....... .. ....... ... ....... . . . ... . . . ... ... ....... . . . . . ... .. ....... . . ... . . . . .. .. ....... . . . . . . ..........................................................................................................................................................................................................................................................................................................................................................................

~n3

ℓ

L

L

ℓ

ℓ

1 v

~n1 -

tO

v

-

~n2

3

L2

Figure 1: LISA configuration and TDI conventions used. In the following we assume a stationary LISA array (cf. Fig. 1), and denote p~i the vectors from the guiding center O to i. We assume here and in the following that the LISA geometry is measured in light-travel time, e.g. L ≡ ˜ ˜ is measured in units of length. L/c, where L GW The single-arm Doppler response to a GW yslr (t) ∼ ∆ν/ν of the light path s → r along arm l can be shown [5, 2] to be h i GW b b b bl ) Ψl (t − k · p~s − Ll ) − Ψl (t − k · p~r ) , (2.2.1) yslr = (1 + slr k · n 6

where

← → bl · h · n 1 n bl . Ψl ≡ 2 1 − (b k·n bl )2

(2.2.2)

Note the geometrical identity p~r − p~s = Ll n bl slr . We further introduce the time-delay operator yslr,d1 ...dm (t) = yslr (t − Ld1 − ... − Ldm ) .

(2.2.3)

Various TDI-observables can be constructed from the single-arm building blocks (2.2.1) with suitable time-delays (2.2.3) as to cancel the (otherwise dominating) laser-noise. One such set of 3 laser-noise free observables is X, Y, Z, defined as X ≡ y132,322 − y123,233 + y231,22 − y321,33 + y123,2 − y132,3 + y321 − y231 , (2.2.4) and Y and Z are given by cyclic permutations of {1, 2, 3}. Geometrically these observables correspond to a ’double-arm’ interferometer, e.g. for X one light-path is 1 → 3 → 1 → 2 → 1 and the second ’arm’ is 1 → 2 → 1 → 3 → 1. To simplify matters, we will here restrict ourselves to work in the longwavelength limit (LWL). The characteristic timescale on which a GW of frequency f is changing is given by h˙ ∼ 2πf O(h), so the characteristic length-scale is λ/2π, the so-called reduced wavelength. The LWL is there˜ λ/2πc, which whould be fore characterized by assuming |~pl | ∼ L = L/c valid for GW frequencies f 1/(2πL) ∼ 10 mHz, assuming an arm-length of L ∼ 5 × 106 km/c ∼ 17 s. We can therefore Taylor-expand in ε ≡ 2πf L 1. The LWL of the single-arm GW responses (2.2.1) is found as GW yslr =−

← → Ll ˙ n bl · h · n bl + O(ε2 ) . 2

(2.2.5)

Note that the Doppler-readouts contain no zero-order contributions in L, i.e. GW yslr = O(ε). In order to express (2.2.4) in the LWL, we also need to expand the time-delays (2.2.3), namely yslr,d1 ...dm (t) = yslr (t) − y˙ slr (t) (Ld1 + ... + Ldm ) + O(ε2 ) .

(2.2.6)

Using this and the symmetry of the first-order term (2.2.5), we can expand X, defined in (2.2.4), in the form (1)

(1)

X = 4L3 y˙ 123 − 4L2 y˙ 231 + O(ε3 ) . 7

(2.2.7)

Plugging in the expansion (2.2.5) of the Doppler readouts, we find the lowestorder contribution as ← → ¨ X (2) = −2L2 L3 (b n2 ⊗ n b2 − n b3 ⊗ n b3 ) : h , (2.2.8) which corresponds to the standard LWL expression for the measured strain h(t) of ground-based detectors, up to a constant pre-factor of −4L2 L3 and ← → ¨ 23 (t), the second time-derivative of h , i.e. we could write X (2) = −4L2 L3 h ← → ← → ← → where h23 (t) ≡ d23 : h and where we defined the detector-tensor dlj as ← → 1 nl ⊗ n bl − n bj ⊗ n bj ) . d lj ≡ (b 2

(2.2.9)

The remaining observables Y, Z are obtained by cyclic index-permutation ← → ← → ¨ X (2) = −4L2 L3 d23 : h , ← → ← → ¨ Y (2) = −4L3 L1 d31 : h , ← → ← → ¨ Z (2) = −4L1 L2 d12 : h ,

(2.2.10a) (2.2.10b) (2.2.10c)

which gives us the explicit relation (to lowest order) between the TDI← → observables X, Y, Z given in the MLDC, and the GW tensor h as used in the F-statistic analysis, as discussed in Sec. 1.1.

2.3

Wide-parameter search grid

For simplicity we used a “foliated” template grid Freq x Sky in the Doppler parameter space ∆λ = {f, α, δ}, consisting of a isotropic sky-grid with stepsizes at the equator: √ 2m dα(0) = dδ = , (2.3.1) (Rorb /c)2π f while for different latitudes we’ll use dα(δ) = dα(0)/ cos(δ), in order to obtain an isotropic sky-grid. The frequency step-size is given by √ 12m df = , (2.3.2) πT where m is the desired maximal mismatch, f is the search-frequency and T the length of observation. The expression for the frequency-resolution is the standard metric frequency stepsize, while the sky-resolution is approximately valid for observation times T & 1/2 year, and can be derived from the orbital phase-metric. 8

2.4

MLDC2 WDB Pipeline

We used an improved version of our pipeline, which is based on the general approach of our pipeline from MLDC1, but which makes more extensive use of the parameter-space metric for finding local maxima and deciding coincidences. Furthermore, this pipeline and its parameter-tuning were specifically designed to reduce secondary maxima and allow mostly primary candidates to pass. This is based on the empirical observation that the location of secondary maxima seems to vary more strongly between ’detectors’ X, Y and Z than the maxima corresponding to primary candidates, especially if an integration time of only T = 1 y is used. The parameters that were used with this pipeline are: grid-mismatch in 1st-stage wide-parameter search m1 = 0.2 1st-stage detection threshold 2Fth = 20 metric sphere for local-maxima mLM = 12.0 coarse-coincidence mismatch mcoinc1 = 0.8 tight-coincidence mismatch mcoinc2 = 0.20 resolution-increase in each zoom zoomFactor = 10 number of final zoom-steps zoomLevel = 2 Using these parameters, the pipeline performs reasonably well in reducing “false alarms” by secondary maxima. Running a search in the frequencyrange f ∈ [5 mHz, 10 mHz] on the training-data set of Challenge-2.1, the pipeline recovers 744 ’correct’ candidates (that lie within m ≤ 0.3 of a signal given in the key) out of the total of 2315 signals in this frequency-range, and produces 7 “false-alarms”, i.e. secondary maxima that managed to cross the pipeline. We therefore expect a false-alarm rate of secondary maxima of roughly ∼ 1% of the number of recovered signals.

References [1] B. Abbott et al. Coherent searches for periodic gravitational waves from unknown isolated sources and Scorpius X-1: results from the second LIGO science run. submitted, 2006. (LIGO Scientific Collaboration) (preprint gr-qc/0605028). [2] J. W. Armstrong, F. B. Estabrook, and M. Tinto. Time-Delay Interferometry for Space-based Gravitational Wave Searches. ApJ, 527:814–826, Dec. 1999. doi: 10.1086/308110.

9

Single−IFO:

TDI X

Single−IFO: TDI Y

Single−IFO:

TDI Z

Wide−parameter F−stat using only Tobs = 1 year



Keep N loudest candidates above threshold 2F > 2Fth

Keep N loudest candidates

Keep N loudest candidates

above threshold 2F > 2Fth

above threshold 2F > 2Fth

Find local maxima within



metric spheres of mismatch



m < mLM

m < mLM

m < mLM

Keep only maxima coincident within metric sphere mCOINC1

ZOOM in coincidences by zoomFactor

Keep only candidates coincident within metric sphere mCOINC2

ZOOM−IN followup by zoomFactor using Tobs = 2 years and coherent

X + (Y−Z)

Repeat zoomLevel times

Figure 2: Schematic representation of wide-parameter pipeline used in MLDC2.

10

[3] K. A. Arnaud, G. Auger, S. Babak, J. Baker, M. J. Benacquista, E. Bloomer, D. Brown, J. B. Camp, J. K. Cannizzo, N. Christensen, J. Clark, N. J. Cornsih, J. Crowder, C. Cutler, S. Fairhurst, L. S. Finn, H. Halloin, K. Hayama, M. Hendry, O. Jeannin, A. Krolak, S. L. Larson, I. Mandel, C. Messenger, R. Meyer, S. Mohanty, R. Nayak, K. Numata, A. Petiteau, M. Pitkin, E. Plagnol, E. K. Porter, R. Prix, C. Roever, B. S. Sathyaprakash, A. Stroeer, P. Sutton, R. Thirumalainambi, D. E. Thompson, J. Toher, R. Umstaetter, M. Vallisneri, A. Vecchio, J. Vinet, J. T. Whelan, and G. Woan. Report on the first round of the Mock LISA Data Challenges. (preprint gr-qc/0701139), 2007. [4] C. Cutler and B. F. Schutz. The generalized F-statistic: multiple detectors and multiple gravitational wave pulsars. Phys. Rev. D., 72:063006, 2005. [5] F. B. Estabrook and H. D. Wahlquist. Response of Doppler spacecraft tracking to gravitational radiation. GReGr, 6:439–447, 1975. [6] P. Jaranowski, A. Królak, and B. F. Schutz. Data analysis of gravitational-wave signals from spinning neutron stars: The signal and its detection. Phys. Rev. D., 58:063001, 1998. [7] A. Królak, M. Tinto, and M. Vallisneri. Optimal filtering of the lisa data. Phys. Rev. D., 70:022003, 2004. [8] LIGO Scientific Collaboration. LAL/LALApps: FreeSoftware (GPL) tools for data-analysis. http://www.lsc-group.phys.uwm.edu/daswg/. [9] R. Prix and J. T. Whelan. F-statistic search for white-dwarf binaries in the first Mock LISA Data Challenge. submitted to CQG, 2007.

11