9512.net

甜梦文库

甜梦文库

当前位置：首页 >> >> # Gravitational waves from coalescing binaries detection strategies and Monte Carlo estimatio

Gravitational waves from coalescing binaries: detection strategies and Monte Carlo estimation of parameters

R. Balasubramanian, B.S. Sathyaprakash, S. V. Dhurandhar

Inter-University Centre for Astronomy and Astrophysics, Post Bag 4, Ganeshkhind, Pune 411 007, India (February 7, 2008) The detection of gravitational waves from astrophysical sources is probably one of the most keenly awaited events in the history of astrophysics. The paucity of gravitational wave sources and the relative di?culty in detecting such waves, as compared to those in the electromagnetic domain, necessitate the development of optimal data analysis techniques to detect the signal, as well as to extract the maximum possible information from the detected signals. Coalescing binary systems are one of the most promising sources of gravitational waves. This is due to the fact that such sources are easier to model and thus one can design detection strategies particularly tuned to such signals. A lot of attention has been devoted in the literature studying such techniques and most of the work has revolved around the Weiner ?ltering and the maximum likelihood estimators of the parameters of the binary system. We investigate such techniques with the aid of di?erential geometry which provides geometric insight into the problem. Such a formalism allows to explore the merits and demerits of a detection scheme independent of the parameters chosen to represent the waveform. The formalism also generalises the problem of choosing an optimal set of templates to detect a known waveform buried in noisy data. We stress the need for ?nding a set of convenient parameters for the waveform and show that even after the inclusion of the second-order post-Newtonian corrections, the waveform can essentially be detected by employing a one-dimensional lattice of templates. This would be very useful both for the purpose of carrying out the simulations as well as for the actual detection process. After setting up such a formalism we carry out a Monte Carlo simulation of the detection process for the initial LIGO/VIRGO con?guration for the ?rst post-Newtonian corrected coalescing binary waveform. We compare the results of our simulations with the currently available estimates of the accuracies in the determination of the parameters and the probability distribution of the maximum likelihood estimators. Our results suggest that the covariance matrix underestimates, by over a factor of two, the actual errors in the estimation of parameters even when the signal-to-noise ratio is as high as 20. As only a tiny fraction of the events is expected to be detected with a signal-to-noise higher than this value, the covariance matrix is grossly inadequate to describe the errors in the measurement of the parameters of the waveform. It is found from our Monte Carlo simulations that the deviations from the covariance matrix are more in the case of the ?rst post-Newtonian waveform than in the case of the Newtonian one. Inclusion of higher-order post-Newtonian corrections introduces new parameters that are correlated with those at the lower post-Newtonian waveform. Such correlations are expected to further increase the discrepancy of the covariance matrix results with those inferred from Monte Carlo simulations. Consequently, numerical simulations that take into account postNewtonian corrections beyond the ?rst post-Newtonian order are needed in order to get a clearer picture about the accuracy in the determination of parameters. We ?nd that with the aid of the instant of coalescence the direction to the source can determined more accurately than with the time of arrival.

arXiv:gr-qc/9508011v1 5 Aug 1995

I. INTRODUCTION

Laser interferometric detectors of gravitational waves such as the LIGO [1] and VIRGO [2] are expected to be operational by the turn of the century. Gravitational waves from coalescing binary systems of black holes and neutron stars are relatively ‘clean’ waveforms in the sense that they are easier to model and for this reason they are amongst the most important candidate sources for interferometric detectors. Binary systems are also valuable sources of astrophysical information as one can probe the universe up to cosmological distances. For instance, statistical analysis of several binary coalescences enables the estimation of the Hubble constant to an accuracy better than 10% [3,4]. Events that produce high signal-to-noise ratio can be potentially used to observe such non-linear e?ects, such as gravitational wave tails, and to put general relativity into test in the strongly non-linear 1

regime [5]. Due to the weak coupling of gravitational radiation with matter the signal waveform is in general very weak and will not stand above the detector noise. In addition to the on-going e?orts to reduce the noise, and hence increase the sensitivity of the detector, a considerable amount of research activity has gone into the development of e?cient and robust data analysis techniques to extract signals buried in very noisy data. For a recent review on gravitational waves from compact objects and their detection see Thorne [6,7]. Various data analysis schemes have been suggested for the detection of the ‘chirp’ waveform from such systems [8–10]. Among them the technique of Weiner ?ltering is the most promising [10–12]. Brie?y, this technique involves correlating the detector output with a set of templates, each of which is tuned to detect the signal with a particular set of parameters. This requires the signal to be known to a high level of accuracy. The fully general relativistic waveform from a coalescing binary system of stars is as yet unavailable. In the absence of such an exact solution, there have been e?orts to ?nd solutions perturbatively. Most of the work done in this area aims at computing the waveform correct to a high degree of accuracy so that theoretical templates computed based on them will obtain a large signal-to-noise ratio (SNR) when correlated with the detector output if the corresponding signal is present. In general, the number of parameters increases as we incorporate the higher order corrections. It is clear that the number of templates depends upon the number of signal parameters. As a consequence, the computing power for an on-line analysis will be greater for a larger number of parameters. In view of this restriction in computing power it is necessary to choose the templates in an optimal manner. This paper in part deals with this question. Investigations till now have been restricted to either choosing a ?nite subset of the signal space as templates [13,14] or choosing templates from the ‘Newtonian’ or the ‘?rst post-Newtonian’ family of waveforms [15–18]. We generalize this problem by using the language of di?erential geometry. We show that it is unnecessary to restrict oneself to templates that are matched exactly to any particular signal. Di?erential geometry has been used in statistics before (see [19] and references therein) and the standard approach is to treat a set of parameterised probability distributions corresponding to a particular statistical model as a manifold. The parameters of the distribution serve as coordinates on this manifold. In statistical theory one frequently comes across the Fisher information matrix whose inverse gives a lower bound for the errors in the estimation of the parameters of a distribution. The Fisher information matrix turns out to be a very natural metric on the manifold of probability distributions and this metric can be used pro?tably in understanding the properties of a particular statistical model. Here in our paper we treat the set of coalescing binary signals corresponding to various parameters of the binary as a submanifold in the linear space of all detector outputs. We show in this paper that both the above mentioned manifolds are equivalent as far as their metrical properties are concerned. The geometric approach turns out to be useful not only in clarifying various aspects of signal analysis but also helps us to pose the question of optimal detection in a more general setting. Once a signal has been detected we can estimate the parameters of the binary. We assume that the parameters of the signal are the same as those of the template with which the maximum correlation is obtained. The errors involved in such an estimation has been worked out by several authors [5,18,20–27], for the case of ‘high’ SNR and for the Newtonian and post-Newtonian waveforms using a single and a network of detectors. For the case of low SNRs one has to resort to numerical simulations. We have started a project to carry out exhaustive numerical simulations speci?cally designed to compute the errors in the estimation of parameters and covariances among them at various post-Newtonian orders, for circular and eccentric orbits, with and without spin e?ects and for di?erent optical con?gurations of the interferometer. In this paper we report the results for the case of the initial LIGO con?guration taking into account only the ?rst post-Newtonian corrections and assuming circular orbits. Going beyond this needs tremendous amount of computing power which is just becoming available. The rest of the paper is organized as follows. In section II we describe the waveform from a coalescing binary system at various post-Newtonian orders. We introduce, following [28], a set of parameters called ‘chirp times’. These parameters are found to be very convenient when we carry out Monte Carlo simulations. It turns out that the covariance matrix is independent of these parameters and hence it is su?cient to carry out the simulations only 2

for a particular set of parameters. In section III we develop a geometric interpretation of signal analysis. We begin by introducing a metric on the manifold from a scalar product, which comes naturally from the theory of matched ?ltering and then show that this metric is the same as the one used by Amari [19]. Using the geometric approach we address the question of optimal ?lter placement and show that for the purpose of detection it is optimal to choose the templates outside the signal manifold. The covariance matrix of errors and covariances is shown to be the inverse of the metric on the manifold. In section IV we discuss the results of our simulations and compare the numerically obtained values and those suggested by the covariance matrix. We ?nd substantial discrepancies in the predictions of the two methods. It is believed that the coalescing binary waveform shuts o? abruptly at the onset of the plunge orbit. This has a major e?ect on the computations of the covariance matrix as well as on the Monte Carlo simulations. We discuss the e?ects of higher post-Newtonian corrections to the waveform. We also emphasisize the use of the instant of coalescence as a parameter in order to determine the direction to the source rather than the time of arrival [29]. Finally in section V we summarise our results and indicate future directions.

II. COALESCING BINARY WAVEFORMS

For the purpose of constructing templates for on-line detection, it is su?cient to work with the so called restricted post-Newtonian gravitational waveform. In this approximation the post-Newtonian corrections are incorporated only in the phase of the waveform, while ignoring corresponding corrections to the amplitude [30]. Consequently, the restricted postNewtonian waveforms only contain the dominant frequency equal to twice the orbital frequency of the binary computed up to the relevant order. In the restricted post-Newtonian approximation the gravitational waves from a binary system of stars, modeled as point masses orbiting about each other in a circular orbit, induce a strain h(t) at the detector given by h(t) = A(πf (t))2/3 cos [?(t)] , (2.1)

where f (t) is the instantaneous gravitational wave frequency, the constant A involves the distance to the binary, its reduced and total mass, and the antenna pattern of the detector [10], and the phase of the waveform ?(t) contains several pieces corresponding to di?erent post-Newtonian contributions which can be schematically written as ?(t) = ?0 (t) + ?1 (t) + ?1.5 (t) + . . . . (2.2)

Here ?0 (t) is the dominant Newtonian part of the phase and ?n represents the nth order post-Newtonian correction to it. In the quadrupole approximation we include only the Newtonian part of the phase given by [10] ?(t) = ?0 (t) = 16πfa τ0 1? 5 f fa

?5/3

+ Φ,

(2.3)

where f (t) is the instantaneous Newtonian gravitational wave frequency given implicitly by t ? ta = τ0 1 ? f fa

?8/3

,

(2.4)

τ0 is a constant having dimensions of time given by τ0 = 5 M?5/3 (πfa )?8/3 , 256 (2.5)

and fa and Φ are the instantaneous gravitational wave frequency and the phase of the signal, respectively, at t = ta . The time elapsed starting from an epoch when the gravitational wave 3

frequency is fa till the epoch when it becomes in?nite will be referred to as the chirp time of the signal. In the quadrupole approximation τ0 is the chirp time. The Newtonian part of the phase is characterised by three parameters: (i) the time of arrival ta when the signal ?rst becomes visible in the detector, (ii) the phase Φ of the signal at the time of arrival and (iii) the chirp mass M = (?3 M 2 )1/5 , where ? and M are the reduced and the total mass of the binary, respectively. At this level of approximation two coalescing binary signals of the same chirp mass but of di?erent sets of individual masses would be degenerate and thus exhibit exactly the same time evolution. This degeneracy is removed when post-Newtonian corrections are included. When post-Newtonian corrections are included the parameter space of waveforms acquires an extra dimension. In this paper we show that even when post-Newtonian corrections up to relative order c?4 , where c is the velocity of light, are included in the phase of the waveform it is possible to make a judicious choice of the parameters so that the parameter space essentially remains only three dimensional as far as the detection problem is concerned. It should, however, be noted that the evolution of the waveform must be known to a reasonably high degree of accuracy and that further o?-line analysis would be necessary to extract useful astrophysical information. With the inclusion of corrections up to second post-Newtonian order the phase of the waveform becomes [31] ?(t) = ?0 (t) + ?1 (t) + ?1.5 (t) + ?2 (t) (2.6)

where ?0 (t) is given by (2.3) and the various post-Newtonian contributions are given by ?1 (t) = 4πfa τ1 1 ? f fa

?1

(2.7)

?1.5 (t) = ?5πfa τ1.5 1 ? and ?2 (t) = 8πfa τ2 1 ? f fa

f fa

?2/3

(2.8)

?1/3

.

(2.9)

Now f (t) is the instantaneous gravitational wave frequency correct up to second postNewtonian order given implicitly by t ? ta = τ0 1 ? f fa

?8/3

+ τ1 1 ?

f fa

?2

? τ1.5 1 ?

f fa

?5/3

+ τ2 1 ?

f fa

?4/3

.

(2.10) In the above equations τ ’s are constants having dimensions of time which depend only on the two masses of the stars and the lower frequency cuto? of the detector fa . The total chirp time now consists of four pieces: the Newtonian contribution τ0 is given by (2.5) and the various post-Newtonian contributions are τ1 = 5 192?(πfa )2 1 8? 743 11 + η , 336 4 m 5 π 2 fa

1/3

(2.11)

τ1.5 = and

(2.12)

4

τ2 =

5 128?

m 2f 2 π a

2/3

3058673 5429 617 2 + η+ η 1016064 1008 144

(2.13)

where η = ?/M. The phase (2.6) contains the reduced mass ? in addition to the chirp mass M. Taking (M, η) to be the post-Newtonian mass parameters the total mass and the reduced mass are given by M = Mη ?3/5 , ? = Mη 2/5 . Note that in the total chirp time τ of the signal the 1.5 post-Newtonian contribution appears with a negative sign thus shortening the epoch of coalescence. With the inclusion of higher order post-Newtonian corrections a chirp template is characterised by a set of four parameters which we shall collectively denote by λ? , ? = 1, . . . , 4. At the ?rst post-Newtonian approximation instead of working with the parameters λ? = {ta , Φ, M, η} we can equivalently employ the set {ta , Φ, τ0 , τ1 } for the purpose of constructing templates. This, as we shall see later, has some advantageous. However, at post-Newtonian orders beyond the ?rst we do not have a unique set of chirp times to work with. The parameters ta and Φ are kinematical that ?x the origin of the measurement of time and phase, respectively, while the Newtonian and the post-Newtonian chirp times are dynamical parameters in the sense that they dictate the evolution of the phase and the amplitude of the signal. It may be mentioned at this stage that in most of the literature on this subject authors use the set of parameters {tC , ΦC , M, η} where tC is the instant of coalescence and ΦC is the phase of the signal at the instant of coalescence. In terms of the chirp times we have introduced, tC is the sum of the total chirp time and the time of arrival and ΦC is a combination of the various chirp times and Φ : tC = ta + τ0 + τ1 ? τ1.5 + τ2 ; ΦC = Φ + 16πfa τ0 + 4πfa τ1 ? 5πfa τ1.5 + 8πfa τ2 . 5 (2.14) (2.15)

In the stationary phase approximation the Fourier transform of the restricted ?rst-postNewtonian chirp waveform for positive frequencies is given by [10,13,21,22]

6

? ) = N f ?7/6 exp i h(f where N = Aπ 2/3

?=1

ψ? (f )λ? ? i

π 4

(2.16)

2τ0 3

1/2 4/3 fa

is a normalisation constant, λ? , ? = 1, . . . , 6, represent the various post-Newtonian parameters λ? = {ta , Φ, τ0 , τ1 , τ1.5 , τ2 } , and ψ1 = 2πf, ψ2 = ?1, ψ3 = 2πf ? (2.18) (2.19) 16πfa 6πfa + 5 5 f fa f fa f fa f fa

?5/3

(2.17)

, ,

?2/3

(2.20) (2.21)

?1

ψ4 = 2πf ? 4πfa + 2πfa ψ5 = ?2πf + 5πfa ? 3πfa ψ6 = 2πf ? 8πfa + 6πfa

,

?1/3

(2.22) (2.23)

.

5

? For f < 0 the Fourier transform is computed using the identity ? h(?f ) = h? (f ) obeyed by real functions h(t). In addition to the above mentioned parameters we shall introduce an amplitude parameter A in section III.

III. A GEOMETRIC APPROACH TO SIGNAL ANALYSIS

In this section we apply the techniques of di?erential geometry to the problem of detecting weak signals embedded in noise. In section III A we introduce the concept of the signal manifold and elaborate on the relationship of our approach with that of Amari [19]. Our discussion of the vector space of all detector outputs is modeled after the discussion given in [32]. In section III B we deal with the problem of choosing a set of ?lters for on-line analysis which would optimise the task of detection of the signal. In section III C we deal with the dimensionality of the chirp manifold when we incorporate higher order post-Newtonian corrections. It is found, that due to covariances between the parameters, it is possible to introduce an e?ective dimension which is less than the dimension of the manifold. This has very important implications for the detection problem.

A. Signal manifold

The output of a gravitational wave detector such as the LIGO, will comprise of data segments, each of duration T seconds, uniformly sampled with a sampling interval of ?, giving the number of samples in a single data train to be N = T /?. Each data train can be considered as a N -tuple (x0 , x1 , . . . , xN ?1 ) xi being the value of the output of the detector at time i?. The set of all such N tuples constitutes a N -dimensional vector space V where the addition of two vectors is accomplished by the addition of corresponding time samples. For later convenience we allow each sample to take complex values. A natural basis for i this vector space is the time basis ei = δm where m and i are the vector and component m indices respectively. Another basis which we shall use extensively is the Fourier basis which ? is related to the time basis by a unitary transformation U : 1 ? ?m = U mn en = √ e N 1 ? em = U ?mn?n = √ e N

N ?1

en exp

n=0 N ?1 n=0

2πimn , N 2πimn . N

(3.1) (3.2)

?n exp ? e

All vectors in V are shown in boldface, and the Fourier basis vectors and components of vectors in the Fourier basis are highlighted with a ‘tilde’. In the continuum case each data train can be expanded in a Fourier series and will contain a ?nite number of terms in the expansion, as the output will be band limited. The expansion is carried out over the exponential functions exp (2πimt/T ) which are precisely the Fourier basis vectors de?ned above. Though the index m takes both positive and negative values corresponding to positive and negative frequencies it is both, possible and convenient to allow m to take only positive values [33]. Thus the vector space V can be considered as being spanned by the N Fourier basis vectors implying immediately that the number of independent vectors in the time basis to be also N . This is the content of the Nyquist theorem which states that it is su?cient to sample the data at a frequency which is twice as large as the bandwidth of a real valued signal, where the bandwidth refers to the range of positive frequencies over which the signal spectrum is non zero. This factor of two does not appear in the vector space picture as we allow in general for complex values for the components in the time basis. A gravitational wave signal from a coalescing binary system can be characterised by a set of parameters λ = (λ0 , λ1 , . . . , λp?1 ) belonging to some open set of the p-dimensional real space Rp . The set of such signals s(t; λ) constitutes a p-dimensional manifold S which is embedded in the vector space V. The parameters of the binary act as coordinates on the 6

manifold. The basic problem of signal analysis is thus to determine whether the detector output vector x is the sum of a signal vector and a noise vector, x = s + n, or just the noise vector, x = n, and furthermore to identify which particular signal vector, among all possible. One would also like to estimate the errors in such a measurement. In the absence of the signal the output will contain only noise drawn from a stochastic process which can be described by a probability distribution on the vector space V. The covariance matrix of the noise C ij is de?ned as, C ij = ni n?j , (3.3)

where an * denotes complex conjugation and an overbar denotes an average over an ensemble. If the noise is assumed to be stationary and ergodic then there exists a noise correlation function K(t) such that Cij = K(|i ? j|?). In the Fourier basis it can be shown that the components of the noise vector are statistically independent [11] and the covariance matrix in the Fourier basis will contain only diagonal terms whose values will be strictly positive: ? Cii = ni n?i . This implies that the covariance matrix has strictly positive eigenvalues. The ? ? ? diagonal elements of this matrix Cii constitute the discrete representation of the power spectrum of the noise Sn (f ). We now discuss how the concept of matched ?ltering can be used to induce a metric on the signal manifold. The technique of matched ?ltering involves correlating the detector output with a bank of ?lters each of which is tuned to detect the gravitational wave from a binary system with a particular set of parameters. The output of the ?lter, with an impulse response q, is given in the discrete case as c(m) 1 =√ N

N ?1

xn q ?n exp [?2πimn/N ] . ? ?

n=0

(3.4)

The SNR (ρ) at the output is de?ned to be the mean of c(m) divided by the square root of its variance: ρ≡ c(m) c(m) ? c(m)

2 1/2

.

(3.5)

By maximising ρ we can obtain the expression for the matched ?lter q(m) matched to a particular signal s(t; λ? ) as q(m) (λ? ) = ?n sn (λ? ) exp [2πimn/N ] ? , ? Cnn (3.6)

where ρ has been maximised at the mth data point at the output and where ? = 1, 2, . . . , p where p is the number of parameters of the signal. We now introduce a scalar product in V. For any two vectors x and y,

N ?1 N ?1 ?1 Cij xi y j = i,j=0 n=0

x, s =

xn s?n (λ? ) ? ? . ? Cnn

(3.7)

In terms of this scalar product, the output of the matched ?lter q, matched to a signal s(λ? ), can be written as, 1 c(m) (λ? ) = √ x, s . N (3.8)

? As Cii is strictly positive the scalar product de?ned is positive de?nite. The scalar product de?ned above on the vector space V can be used to de?ne a norm on V which in turn can be used to induce a metric on the manifold. The norm of a vector x is de?ned as 1/2 1/2 x = x, x . The ratio for the matched ?lter can be calculated to give ρ = s, s . The √ 1/2 norm of the noise vector will be a random variable n, n with a mean of N as can be 7

seen by writing the expression for the norm of the noise vector and subsequently taking an ensemble average. The distance between two points in?nitesimally separated on S can be expressed as a quadratic form in the di?erences in the values of the parameters at the two points: g?ν dλ? dλν ≡ s(λ? + dλ? ) ? s(λ? ) = = ?s ?s , ?λ? ?λν dλ? dλν ?s dλ? ?λ? (3.9) (3.10)

The components of the metric in the coordinate basis are seen to be the scalar products of the coordinate basis vectors of the manifold. Since the number of correlations we can perform on-line is ?nite, we cannot have a ?lter corresponding to every signal. A single ?lter though matched to a particular signal will also detect signals in a small neighbourhood of that signal but with a loss in the SNR. The metric on the manifold quanti?es the drop in the correlation in a neighbourhood of the signal chosen. Taking the output vector to be x and two signal vectors s(λ) and s(λ + dλ) and using Schwarz’s inequality we have, x, s(λ + dλ) ? x, s(λ) = x, s(λ + dλ) ? s(λ) ≤ x s(λ + dλ) ? s(λ) = x g?ν dλ? dλν . (3.11) (3.12)

As is apparent the drop in the correlation can be related to the metric distance on the manifold between the two signal vectors. We now discuss Amari’s [19] work in the context of using di?erential geometry in statistics and elaborate on the relationship with the approach we have taken. The set of parametrised probability distributions corresponding to a statistical model constitute a manifold. The parameterized probability distributions in the context of signal analysis of gravitational waves from coalescing binaries are the ones which specify the probability that the output vector will lie in a certain region of the vector space V given that a signal s(t; λ) exists in the output which we denote as p(x|s(t; λ)) Since it is not our intention to develop Amari’s approach any further we will be brief and will make all the mathematical assumptions such as in?nite di?erentiability of functions, interchangeability of the di?erentiation and expectation value operators, etc. The set of probability distributions p(x|s(λ)) where λ ∈ Rp , constitutes a manifold P of dimension p. At every point on this manifold we can construct a tangent space T 0 on ? which we can de?ne the coordinate basis vectors as ?? = ?λ? . Any vector A in this tangent space can be written as a linear combination of these coordinate basis vectors. We now ? de?ne p random variables σ? = ? ?λ? log(p(x|s(λ))). It can easily be shown that σ? = 0. We assume that these p random variables are linearly independent. By taking all possible linear combinations of these random variables we can construct another linear space T 1 . Each vector B in T 1 can be written as B = B ? σ? . The two vector spaces T 0 and T 1 are isomorphic to each other, which can be shown explicitly by making the correspondence σ? ? ?? . The vector space T 1 has a natural inner product de?ned on it which is the covariance matrix of the p random variables σ? . This scalar product can be carried over to T 0 using the correspondence stated above. The metric on the manifold can be de?ned by taking the scalar product of the coordinate basis vectors g?ν = ?? , ?ν = σ? σν (3.13)

In statistical theory the above matrix g?ν is called the Fisher information matrix. We will also denote the Fisher matrix, as is conventional, as Γ?ν . It is clearly seen that orthogonality between vectors in the tangent space of the manifold is related to statistical independence of random variables in T 1 . If we take the case of Gaussian noise the metric de?ned above is identical to the one obtained on the signal manifold by matched ?ltering. Gaussian noise can be described by the distribution,

8

p(n) =

1 exp ? 2

N ?1 j,k=0

?1 Cjk nj nk? 1/2

[ (2π)N det [Cij ] ]

=

1 exp ? 2

N ?1 j,k=0

? ?1 ? ? Cjk nj nk?

1/2

=

1 exp ? 2

N ?1 i=0

ni ni? ? ? ? Cii 1/2

,

? (2π)N det Cij

? (2π)N det Cij

(3.14) ? where in the last step we have used the diagonal property of the matrix C ij which implies ? ?1 = 1/Cii . ? that Cii As the noise is additive p(x|s(λ)) can be written as p(x ? s(λ)). Assuming Gaussian noise we can write the expressions for the random variables σi as, σ? = 1 2 ? x ? s(λ), x ? s(λ) = ?λ? n, ? s(λ) , ?λ? (3.15)

where in the last step we have used x = s(λ) + n. The covariance matrix for the random variables σ? can be calculated to give σ? σν = ? ? s(λ), ν s(λ) , ?λ? ?λ (3.16)

which is the same metric as de?ned over the signal manifold. Thus, both the manifolds S and P are identical with respect to their metrical properties. We will henceforth restrict our attention to the signal manifold S. For the purpose of our analysis we will choose a minimal set of parameters characterizing the gravitational wave signal from a coalescing binary. We consider only the ?rst post-Newtonian corrections. In section II we have already introduced the four parameters λ? = {ta , Φ, τ0 , τ1 }. We now introduce an additional parameter for the amplitude and call it λ0 = A. The signal can now be written as s(f ; λ) = Ah(f ; ta , Φ, τ0 , τ1 ), where, ? λ ≡ {A, ta , Φ, τ0 , τ1 }. Numerically the value of the parameter A will be the same as that of the SNR obtained for the matched ?lter. We can decompose the signal manifold into a manifold containing normalised chirp waveforms and a one-dimensional manifold corresponding to the parameter A. The normalised chirp manifold can therefore be parameterized by ? {ta , Φ, τ0 , τ1 }. This parameterization is useful as the coordinate basis vector ?A will be orthogonal to all the other basis vectors as will be seen below. In order to compute the metric, and equivalently the Fisher information matrix, we use the continuum version of the scalar product as given in [20], except that we use the two sided power spectral density. This has the advantage of showing clearly the range of integration in the frequency space though we get the same result using the discrete version of the scalar product. Using the de?nition of the scalar product we get

∞

g?ν =

fa

s s df ??(f ; λ) ??? (f ; λ) + c.c. ? Sn (f ) ?λ ?λν

(3.17)

Recall that in the stationary phase approximation the Fourier transform of the coalescing ? ? ? binary waveform is given by h(f ) = N f ?7/6 exp i ? ψ? (f )λ? and s(f ) = Ah(f ), where ψ? (f ) are given by equations (2.18-2.23), ? = 1, . . . , 4 and λ? = {ta , Φ, τ0 , τ1 }. Note, in particular, that in the phase of the waveform the parameters occur linearly thus enabling a very concise expression for the components of g?ν . The various partial derivatives are given by ??(f ; λ) s = iψ? (f )?(f ; λ), s ?λ? (3.18)

where we have introduced ψ0 = ?i/A. On substituting the above expressions for the partial derivatives in eq. (3.17) we get,

∞ ? ? ψ? (f )ψν (f ) h(f ) 2

g?ν = ψ? h, ψν h = 2

fa

Sn (f ) 9

df

(3.19)

The above de?nition of the amplitude parameter A, as in Culter and Flannagan [22], disjoins the amplitude of the waveform from the rest of the parameters. Since ψ0 is pure imaginary and ψ? ’s are real, it is straightforward to see from eq. (3.19) that g00 = 1 and g0? = 0; ? = 1, . . . , 4. (3.20)

The rest of the components g?ν are seen to be independent of all the parameters except A i.e. g?ν ∝ A2 . As A is unity for the normalised manifold the metric on the normalised manifold is ?at. This implies not only that the manifold is intrinsically ?at (in the stationary phase approximation) but also that the coordinate system used is Cartesian. If instead of the chirp time τ0 we use the parameter M then the metric coe?cients will involve that parameter and the coordinate system will no longer remain Cartesian.

B. Choice of ?lters

We now use the formalism developed to tackle the issue of optimal ?lter placement. Till now, it has been thought necessary to use a ?nite subset of the set of chirp signals as templates for detection. We show that this is unduly restrictive. We suggest a procedure by which the detection process can be made more ‘e?cient’ by moving the ?lters out of the manifold. It must be emphasised that the algorithm presented below is both, simplistic and quite adhoc and is not necessarily the best. Moreover, we have implemented the algorithm only for the Newtonian case where the computational requirements are not very heavy. The signal manifold corresponding to higher post-Newtonian corrections will be a higher dimensional manifold and here the computational requirements will be substantial. The choice of optimal ?lters which span the manifold will then be crucial. Detection of the coalescing binary signal involves computing the scalar product of the output of the detector with the signal vectors. Subsequently one would have to maximise the correlations over the parameters and the number so obtained would serve as the statistic on the basis of which we can decide whether a signal is present in the given data train. Geometrically, this maximization corresponds to minimizing the angle between the output vector and the vectors corresponding to the normalised signal manifold. Using the cosine formula, cos(θ) = s(λ), x s(λ) = x s(λ)

2

+ x 2 ? x ? s(λ) 2 x s(λ)

2

,

(3.21)

it is clear that as s(λ) is unity for the vectors belonging to the normalised signal manifold, maximising the scalar product is equivalent to minimising x ? s(λ) which is the distance between the tip of the output vector and the manifold. Given the constraints of computational power one would be able to evaluate only a ?nite number of these scalar products, say nF , in a certain amount of time depending on the data train length and the padding factor. It is therefore necessary to be able to choose the nF ?lters in such a manner that the detection probability is maximal. We will need e?cient on-line data analysis for two reasons: (i) To isolate those data trains which have a high probability of containing a signal and (ii) to determine the parameters of the binary early on during the inspiral and to use them for dynamical recycling techniques. Due to the ?niteness of the ?lter spacing the signal parameters will in general not correspond to any of the nF ?lters chosen and this will lead to a drop in the maximum possible correlation. Till now attention has been focussed on identifying a set of optimal set of ?lters which are a discrete subset of the manifold. If detection is the sole purpose then the di?erential geometric picture suggests that con?ning the ?lter vectors to the signal manifold is an unnecessary restriction and in fact non optimal. Thus it is worthwhile to explore making a choice of ?lters outside the manifold. The ?lter vectors will thus belong to V but will not, in general, correspond to any signal. It is, of course, true that we are sacri?cing on the maximum possible correlation obtainable (when the signal’s parameters coincide with those of the ?lter). Thus the problem essentially is to select nF ?lter vectors which optimize the detection the e?ciency of which depends upon the properties of the manifold. 10

In general a single ?lter vector would have to pick up signals over a region of the manifold. The extent of this region is determined by ?xing a threshold on the correlation between the ?lter and any signal in the region. We will denote this threshold by η, where η takes a value which is close to, but less than unity. The typical value suggested for η is ? 0.8, [13]. For a given ?lter q and a threshold η the region on the manifold corresponding to the ?lter will be denoted as Sq (η), where Sq (η) ? S. Geometrically this region is the intersection of an open ball in V of radius 21/2 (1 ? η)1/2 (using the distance de?ned by the scalar product), with center q, and the manifold. The nF ?lters taken together would have to ‘span’ the manifold which means that the union of the regions covered by each ?lter would be the manifold itself i.e., q Sq (η) = S. If the ?lter q lies on the manifold, then the correlation function cq (λ) = s(λ), q will reach its maximum value of unity in Sq (η) when q = s(λ) and will fall o? in all directions. This means that the signals in the region which are further away from the ?lter are less likely to be picked up as compared to those in the immediate neighbourhood of the ?lter q. We assume that a ?nite subset of the normalised signal manifold has been chosen to act as ?lters by some suitable algorithm [13], which taken together span the manifold. The number of ?lters will be determined by the available computing power. Consider one of these ?lters q, the region corresponding to it for a threshold of η, Sq (η), and an arbitrary normalised vector qo which belongs to V but not necessarily S. By correlating the vector qo with vectors in Sq (η) we obtain the correlation function cqo (λ) = s(λ), qo . We select that qo to serve as a more optimal ?lter which maximises the average of this correlation function: s(λ), qo

av

=

Sq (η)

1 √ p gd λ

s(λ), qo

Sq (η)

√ p gd λ =

Sq (η)

s(λ)dp λ , qo

,

(3.22)

where g = det [g?ν ]. In the last step above the integration and the scalar product operations in the above equation have been interchanged. Moreover, for the normalised chirp manifold the metric does not depend upon the parameters in the coordinate system we have chosen √ and therefore, g?ν is a constant and the factor g cancels. We now use Schwarz’s inequality to maximise the average correlation to obtain, qo = N

Sq (η)

s(λ)dp λ,

(3.23)

where N is just a normalisation constant. We implemented the above algorithm for ?lter placement for the case of Newtonian signals with certain modi?cations. The normalised chirp waveform consists of three parameters (Φ, ta , τn ). If we keep ta and τn ?xed then the tip of the signal vector traces out a circle as we vary Φ. As any circle lies on a plane we can express a signal vector as a linear sum of two vectors where the two vectors di?er only in the phase parameter and we take this phase di?erence to be π/2. Thus we need only two mutually orthogonal ?lters to span the phase parameter. The time of arrival parameter ta is also a ‘convenient parameter’ as by the use of fast Fourier transforms the correlations for arbitrary time of arrivals can be performed at one go. It is therefore not pro?table for us to maximise the average correlation over the phase parameter Φ and the time of arrival ta . In view of the above restrictions we modi?ed the the ?lter placement algorithm. We consider the correlation function for the case when the ?lter vector is on the manifold. We de?ne the ‘line of curvature’ to be the curve on the manifold along which the correlation function falls the least. Figure (1) illustrates the correlation function plotted as a function of τ0 along the line of curvature. It is seen from the contour diagram of the numerically computed correlation function that the line of curvature lies nearly on the submanifold ta + τn = constant, of the normalised chirp manifold. We take two curves passing through the point q in the region Sq (η): 1. ta + τn = constant, Φ = 0 and 2. ta + τn = constant, Φ = π/2. 11

We obtain one ?lter for each of the two curves by evaluating eq. 3.23 where the domains of integration correspond to the segments of the curves de?ned. Having determined the two ?lters we again plot the correlation function along the line of curvature as a function of τ0 in ?gure (1). The region of the manifold selected corresponds to a range of 5.8 secs to 6.0 secs in the parameter τ0 . It is interesting to note that we could have translated these values keeping the di?erence same without a?ecting our results. It can be seen that the correlation has a minimum at the center. In order to get a ?atter correlation curve we select a linear combination of the original ?lter and the one obtained by averaging with suitable weights attached to each ?lter. This performs reasonably well as shown by the thick curve in the ?gure. The importance of having a ?atter correlation function lies in the fact all the signals in a region can be picked up with equal e?ciency and the drop in the maximum possible correlation can be compensated for by lowering the threshold. The average correlation obtained for the optimal ?lter is only marginally better than that obtained for the ?lter placed on the manifold. In the discussion above, we had started with a ?xed number of ?lters nF on the manifold and obtained another set of nF ?lters which performs marginally better than the former set. Equivalently we can try to increase the span of each ?lter retaining the same threshold but reducing the number of ?lters required. In Figure 1 we observe that the optimal ?lter chosen spans the entire region considered with a threshold greater than 0.9, whereas the ?lter on the manifold spans about half the region at the same threshold. This indicates that by moving the ?lters out of the manifold in the above manner it may be possible to reduce the number of ?lters by a factor of two or so. One must however, bear in mind that the bank of ?lters obtained in this way are not optimal. There is scope to improve the scheme further. Our analysis is indicative of this feature.

C. E?ective dimensionality of the parameter space of a second order post-Newtonian waveform

It has already been shown that the ?rst post-Newtonian waveform is essentially onedimensional [28]. We argue in this subsection that even the second post-Newtonian waveform is essentially one-dimensional and a one-dimensional lattice su?ces to ?lter the waveform. A Newtonian waveform is characterised by a set of three parameters consisting of the time of arrival, the phase of the signal at the time of arrival and chirp mass (equivalently, Newtonian chirp time). In this case, for the purpose of detection, one essentially needs to employ a one-dimensional lattice of ?lters corresponding to the chirp mass, the time of arrival being taken care by the fast Fourier transform algorithm and the phase being determined using a two-dimensional basis of orthogonal templates. When post-Newtonian corrections are included in the phase of the waveform the number of parameters increases apparently implying that one needs to use a two-dimensional lattice of ?lters corresponding to, say, the chirp and reduced masses (equivalently the Newtonian and post-Newtonian chirp times) which in turn means that the number of templates through which the detector output needs to be ?ltered goes up by several orders of magnitude. One of us (BSS) has recently shown that for the purpose of detection it is su?cient to use a one-dimensional lattice of ?lters even after ?rst post-Newtonian corrections are included in the phase of the waveform and the relevant parameter here is the sum of the Newtonian and post-Newtonian chirp times. What happens when corrections beyond the ?rst post-Newtonian order are incorporated in the phase of the waveform? The coalescing binary waveform is now available up to second post-Newtonian order [31,34]. Blanchet et al. argue that the phase correction due to the second order postNewtonian (2PN) term induces an accumulated di?erence of 10 cycles in a total of 16000. Consequently, it is important to incorporate the 2PN terms in the templates. When the 2PN terms are included it is useful to consider that the full waveform is parameterised by three additional parameters, corresponding to the chirp times at the 1, 1.5 and 2PN order (cf. Sec II), as compared to the Newtonian waveform. Of course, as far as the detection problem is concerned there is only one additional parameter since the chirp times are all functions of the two masses of the binary. However, for the purpose of testing general

12

relativity one can consider each of the chirp times to be independent of the rest [5,23]. Our problem now is to ?nd the dimensionality of the parameter space of a 2PN waveform. To ′ this end we consider the ambiguity function C(λ , λ) which is nothing but the correlation function of two normalised waveforms one of whose parameters (λ) are varied by holding ′ the parameters of the other ?xed (λ ) : C(λ , λ) = q(λ ), q(λ) ;

′ ′ ′

q(λ ), q(λ ) = q(λ ), (λ) = 1.

′

′

(3.24)

It is useful to think of λ to be the parameters of a template and λ to be that of a signal. With this interpretation the ambiguity function simply gives the span of a ?lter in the parameter space. The ambiguity function for the full waveform is a four-dimensional surface since there are four independent parameters. To explore the e?ective dimensionality of the parameter we consider the set of parameters to be {ta , Φ, m1 , m2 }, where m1 and m2 are the two masses of the binary. We have shown the contours of the ambiguity function maximised over ta and Φ (since these two parameters do not explicitly need a lattice of templates) in Fig. 2. The template at the centre of the plot corresponds to a binary waveform with m1 = m2 = 1.4M⊙ and the signal parameters are varied over the entire astrophysically interesting range of masses: m1 , m2 ∈ [1.4, 10]M⊙. From these ?gures we ?nd that the ambiguity function is almost a constant along a particular line in the m1 –m2 plane. This means that a template at the centre of the grid spans a relatively large area of the parameter space by obtaining a correlation very close to unity for all signals whose masses lie on the curve along which the ambiguity function roughly remains a constant. It turns out that the equation of this curve is given by τ0 + τ1 ? τ1.5 + τ2 = const. (3.25)

Let us suppose we begin with a two-dimensional lattice of ?lters corresponding to a certain grid (albeit, nonuniform) laid in the m1 –m2 plane. Several templates of this set will have their total chirp time the same. Now with the aid of just one template out of all those that have the same chirp time we can e?ectively span the region that is collectively spanned by all such ?lters. More precisely, we will not have an appreciable loss in the SNR in replacing all templates of a given total chirp time by one of them. Consequently, the signal manifold can be spanned by a one-dimensional lattice of templates.

IV. ESTIMATION OF PARAMETERS

In this section we discuss the accuracy at which the various parameters of a coalescing binary system of stars can be estimated. All our results are for a single interferometer of the initial LIGO-type which has a lower frequency cuto? at 40 Hz. At present it is beyond the computer resources available to us to carry out a simulation for the advanced LIGO. In the ?rst part of this section we brie?y review the well known results obtained for the variances and covariances in the estimation of parameters using analytical methods. Analytical methods assume that the SNR is su?ciently large (the so-called strong signal approximation) and implicitly use a continuum of the parameter space. In reality, however, these assumptions are not necessarily valid and hence it is essential to substantiate the results obtained using analytical means by performing numerical simulations. In the second part of this section we present an exhaustive discussion of the Monte Carlo simulations we have performed to compute the errors and covariances of di?erent parameters. As we shall discuss below the computation of errors using the covariance matrix is erroneous even at an SNR of 10-20. Our estimation of 1-σ uncertainty in the various parameters, at low SNRs, is substantially larger than those computed using the covariance matrix. However, for high values of the SNR (> 25–30) Monte Carlo estimation agrees with the analytical results.

13

A. Covariance matrix

In recent years a number of authors have addressed issues related to the variances expected in the parameter estimation [5,18,20–27], In the standard method of computing the variances in the estimation of parameters one makes the assumption that the SNR is so large that with the aid of such an approximation one can ?rst construct the Fisher information matrix Γ?ν and then take its inverse to obtain the covariance matrix C?ν . In the strong signal approximation the Fisher information matrix and the covariance matrix are given by g?ν = Γ?ν = ?s ?s , ?λ? ?λν ; C?ν = Γ?1 ?ν. (4.1)

As we have seen before the Fisher information and consequently the covariance matrix is block diagonal and hence there is no cross-talk, implying vanishing of the covariances between the amplitude and the other parameters. Consequently, we need not construct, for the purpose of Weiner ?ltering, templates corresponding to di?erent amplitudes. The judicious choice of parameters has also allowed a very elegant expression for the Fisher information matrix. It is particularly interesting to note that the information matrix does not depend on the values of the various parameters, except for the amplitude parameter, and hence is a constant as far as the parameter space of the normalised waveforms is considered. For the purpose of numerical simulations it is convenient to choose the set λ? = {A, ta , Φ, τ0 , τ1 } where A is the amplitude parameter, ta and Φ are the time of arrival of the signal and its phase at the time of arrival, and τ0 and τ1 are the Newtonian and the post-Newtonian coalescence times. For noise in realistic detectors, such as LIGO, the elements of the Fisher information matrix cannot be expressed in a closed form and, for the set of parameters employed, it is not useful to explicitly write down the covariance matrix in terms of the various integrals since the errors and covariances do not have any dependence on the parameters. We thus evaluate the information matrix numerically and then take its inverse to obtain the covariance matrix. Instead of dealing with the covariance matrix C it more instructive to work with the matrix of standard deviations and correlation coe?cients D which is related to the former by D?ν = C?ν , if ? = ν C?ν /(σ? σν ) if ? = ν, (4.2)

where σ? = D?? is the 1-σ uncertainty in the parameter λ? . The o?-diagonal elements of D take on values in the range [?1, 1] indicating how two di?erent parameters are correlated: For ? = ν, D?ν = 1, indicates that the two are perfectly correlated, D?ν = ?1 means that they are perfectly anticorrelated and D?ν = 0 implies that they are uncorrelated. Since the information matrix is block-diagonal, the amplitude parameter is totally uncorrelated with the rest and thus an error in the measurement of A will not re?ect itself as an error in the estimation of the other parameters and vice versa. In contrast, as we shall see below, Newtonian chirp time is strongly anticorrelated to post-Newtonian chirp time, which implies that if in a given experiment τ0 happens to be estimated larger than its true value then it is more likely that τ1 will be estimated to be lower than its actual value. Such correlations are useful as far as detection is concerned since they tend to reduce the number of templates needed in ?ltering a given signal. On the other hand, strong correlations increase the volume in the parameter space to which an event can be associated at a given con?dence level. It seems to be in general true that a given set of parameters do not satisfy the twin properties of having small covariances and reducing the e?ective dimension of the manifold for the purpose of ?ltering. We elaborate on this point below. Given a region in a parameter space, it is useful to know the proper volume (as de?ned by the metric) of the manifold corresponding to the said region. In choosing a discrete set of ?lters for the detection problem one has to decide upon the maximum allowable drop in the correlation due to the ?nite spacing. Once this is ?xed, the number of ?lters can be determined from the total volume of the manifold. For the detection problem it is bene?cial to have a small volume whereas if the waveform is parameterized in such a way such that the manifold corresponding to it covers a large volume, then one can determine the parameters 14

to a greater accuracy. As a simple example let us consider a two-dimensional toy model: λ = {λ1 , λ2 }. We compare di?erent signal manifolds each corresponding to a di?erent parameterizations of the waveform. We assume that the covariance matrix and its inverse, the Fisher information matrix, to be: C?ν = σ11 σ12 σ12 σ22 and Γ?ν = γ11 γ12 γ12 γ22 = 1 2 (σ11 σ22 ? σ12 ) σ22 ?σ12 ?σ12 σ11 . (4.3)

The volume of the manifold corresponding to a region K of the parameter space is given as, VK =

K

γ11 γ22 1 ?

2 γ12 dλ1 dλ2 = γ11 γ22

K

γ11 γ22 [1 ? ?] dλ1 dλ2 ,

(4.4)

2 where ? = γ12 /(γ11 γ22 ) is the correlation coe?cient. It can be clearly seen that if the correlation coe?cient is small, keeping the variances in the parameters, ?xed then the volume of the manifold is maximal. Since the parameters τ0 and τ1 are highly anticorrelated the proper volume corresponding to the region reduces to zero showing that the e?ective dimensionality of the manifold is less. Though, in principle, the variances and covariances are independent of the chirp time, in reality there arises an indirect dependence since one terminates a template at a frequency f = 1/(63/2 πM ) (where M is the total mass of the binary) corresponding to the plunge radius at a = 6M. Therefore, larger mass binaries are tracked over a smaller bandwidth so much so that there is less frequency band to distinguish between two chirps of large, but di?erent, total mass. Consequently, at a given SNR the error in the estimation of chirp times is larger for greater mass binaries. This is re?ected by the fact that the integrals in eq. (3.17) are somewhat sensitive to the value of the upper cuto?. (This also explains why the errors in the estimation of the chirp and reduced masses are larger for greater mass binaries [21,22].) In the following we assume that the noise power spectral density is that corresponding to the initial LIGO for which a ?t has been provided by Finn and Cherno? [21]. For an SNR of 10 the matrix D is given by

? D?ν = ?

?

1.0

? 0 0 0 8.37 0.999 ?0.999 ? , 3.16 ?0.998 ? 8.4

(4.5)

for the ?rst post-Newtonian corrected signal. While computing varianaces and covariances the integrals in equation (3.19) are evaluated by chosing a ?nite upper limit of 1 kHz. In the above matrices the entries are arranged in the order {A, ta , Φ, τ0 } in the Newtonian case, {A, ta , Φ, τ0 , τ1 } in the post-Newtonian case, o? diagonal elements are dimensionless correlation coe?cients and, where appropriate, diagonal elements are in ms. The values quoted in the case of the Newtonian waveform are consistent with those obtained using a di?erent set of parameters by Finn and Cherno? [21]. In Fig. 3 we have plotted σ’s, at an SNR of 10, as a function of the upper frequency cuto? fc for Newtonian and post-Newtonian chirp times and the instant of coalescence tC . We see that σ is larger for higher mass binaries but this is because we have ?xed the SNR. However, if we consider binaries of di?erent total masses, all located at the same distance, then a more massive binary produces a stronger SNR so that in reality it may be possible to determine its parameters more accurately than that of a lighter binary. In Fig. 4 we have plotted σ’s for binaries all located at the same distance as a function of total mass. We ?x one of the masses at a value of 1.4 M⊙ and vary the other from 1.4 M⊙ to 10 M⊙ . In computing the σ’s plotted in this ?gure we have 15

for the Newtonian signal, and by ? ? 1.0 0 0 0 0 20.4 0.997 ?0.972 0.911 ? ? ? ? 6.7 ?0.954 0.881 ? , D?ν = ? ? 45.1 ?0.982 ? 25.98

(4.6)

terminated the waveform at the plunge orbit and normalised the SNR of a 10M⊙ -1.4M⊙ binary system to 10. As a function of M the uncertainties in τ0 and τ1 initially fallo? since the increase in the SNR for larger mass binaries more than compensates for the drop in the upper frequency cuto?. However, for M larger than a certain M0 the increase in SNR is not good enough to compensate for the drop in fc so much so that the uncertainties in τ0 and τ1 increase beyond M0 . The parameter ta , however, falls o? monotonically.

B. Monte Carlo estimation of parameters

In this Section we present the ?rst in a series of e?orts to compute the covariance matrix of errors through numerical simulations for a coalescing binary waveform at various postNewtonian orders. Analytical computation of the covariance matrix, as in the previous section, gives us an idea of the covariances and variances but, as we shall see in this Section, at low SNR’s it grossly underestimates the errors. Quite apart from the fact that the assumptions made in deriving the covariance matrix might be invalid at low SNR’s, in a realistic detection and data analysis, other problems, such as discreteness of the lattice of templates, ?nite sampling of the data, etc., do occur. It therefore seems necessary to check the analytical calculations using numerical simulations to gain further insight into the accuracy at which physical parameters can be measured. This Section is divided into several parts: In the ?rst part we highlight di?erent aspects of the simulation, in the second part we brie?y discuss the choice of templates for the simulation, in the third we elaborate on the Monte Carlo method that we have adopted to carry out our simulations and in the fourth we discuss problems that arise in a numerical simulation. The results of our study are discussed in the next Section.

1. Parameters of the simulation

? Let s(t) be a signal of strength A characterised by a set of parameters λ ? ? s(t; λ) = Ah(t; λ), h, h = 1.

k

(4.7)

In data analysis problems one considers a discrete version {s |k = 0, . . . , N ? 1} of the waveform s(t) sampled at uniform intervals in t : sk ≡ s(k?); k = 0, . . . , N ? 1, (4.8)

where ? denotes the constant interval between consecutive samples and N is the total number of samples. The sampled output xk of the detector consists of the samples of the noise plus the signal: xk = nk + sk .

?1

(4.9)

The sampling rate, fs = ? (also referred to as the sampling frequency) is the number of samples per unit time interval. In a data analysis problem the sampling frequency is determined by the signal bandwidth. If B is the signal bandwidth, i.e., if the Fourier transform of the signal is only nonzero over a certain interval B, then it is su?cient to sample at a rate fs = 2B. In our case there is a lower limit in the frequency response of the detector since the detector noise gets very large below a seismic cuto? at about 10–40 Hz. As mentioned in the last Section there is also an upper limit in frequency up to which a chirp signal is tracked since one does not accurately know the waveform beyond the last stable orbit of the binary. This corresponds to gravitational wave frequency fc = 1/(63/2 πM ). For a neutron star-neutron star (ns-ns) binary fc ? 1525 Hz while for a ns-black hole (of 10 M⊙ ) (ns-bh) binary fc ? 375 Hz. Due to constraints arising out of limited computational power, we terminate waveforms at 750 Hz even when fc is larger than 1000 Hz. Such a shuto? is not expected to cause any spurious results since, even in the case of least massive binaries of ns-ns, which we consider in this study, more than 99 % of the ‘energy’ is extracted by the time the signal reaches 750 Hz. We have carried out simulations with two types of upper cuto?: 16

1. one in which all templates, irrespective of their total mass, are shuto? beyond 750 Hz., 2. a second in which the upper frequency cuto? is chosen to be 750 Hz or fc , whichever is lower. Consistent with these cuto?s the sampling rate is always taken to be 2 kHz. (We have carried out simulations with higher sampling rates and found no particular advantage in doing so nor did we ?nd appreciable changes in our results.) In all our simulations, as in the previous section, we take the detector noise power spectral density S to be that corresponding to initial LIGO [21]. For the purpose of simulations we need to generate noise corresponding to such a power spectrum. This is achieved by the following three steps: 1. generate Gaussian white noise n′k with zero mean and unit variance, n′k = 0, n′k n′l = δ kl where an overbar denotes average over an ensemble, 2. compute its Fourier transform 1 n′k ≡ √ ? N and 3. multiply the Fourier components by the square root of the power spectral density, √ ? nk = S k n′k . ? The resultant random process has the requisite power spectrum. In the above, the second step can be eliminated since the Fourier transform of a Gaussian random process is again a Gaussian, but with a di?erent variance. In other words we generate the noise directly in the Fourier domain. The simulated detector output, in the presence of a signal sk , in the Fourier domain is given by xk = nk + sk ? ? ? where sk is the discrete Fourier transform of the signal. ?

2. Choice of templates

N ?1

n′l exp(2πikl/N )

l=0

(4.10)

To ?lter a Newtonian signal we employ the set of parameters {ta , Φ, τ0 } and to ?lter a post-Newtonian signal we employ the set {ta , Φ, τ0 , τ1 }. Templates need not explicitly be constructed for the time of arrival since computation of the scalar product in the Fourier domain (and the availability of fast Fourier transform (FFT) algorithms) takes care of the time of arrival in essentially one computation (N log2 N operations as opposed to N 2 operations, where N is the number of data points). Moreover, there exists a two-dimensional basis for the phase parameter which allows the computation of the best correlation with the aid of just two ?lters. Consequently, the parameter space is essentially one-dimensional in the case of Newtonian signals and two-dimensional in the case of post-Newtonian signals. (However, as shown in Section III C it is to be noted that for the purpose of detection the e?ective dimensionality of the parameter space, even with the inclusion of second post-Newtonian corrections, is only one-dimensional.) We adopt the method described in Sathyaprakash and Dhurandhar [13] to determine the templates needed for chirp times. As described in [13,28] ?lters uniformly spaced in τ0 and τ1 covers the parameters space e?ciently.

17

3. Monte Carlo method

In order to compute variances and covariances numerically, we employ the Monte Carlo method. The basic idea here is to mimic detection and estimation on a computer by performing a very large number of simulations so as to minimize the uncertainties induced by noise ?uctuations. In our simulations we generate a number of detector outputs {xk } each ? corresponding to a de?nite signal sk (λ) of a certain strength, but corresponding to di?erent realisations of the random process {nk }. Computation of the covariance matrix involves ?ltering each of these detector outputs through an a priori chosen set (or lattice) of templates {q(t; t λk )|k = 1, . . . , nf }, where nf denotes the number of templates. The templates of the lattice each has a distinct set of values of the test parameters t λk and together they span a su?ciently large volume in the parameter space. The simulated detector output is correlated with each member of the lattice to obtain the corresponding ?ltered output C(t λk ) : C(t λk ) = x, q(t λk ) . (4.11)

For a given realisation of noise a particular template obtains the largest correlation and its parameters are the measured values m λ of the signal parameters. Thus, the measured values of the parameters are de?ned by max C(t λk ) = C(m λ).

k

(4.12)

The measured values, being speci?c to a particular realisation of noise, are random variables. Their average provides an estimation e λ of the true parameter values and their variance is a measure of the error σλ in the estimation:

eλ

= m λ,

2 σλ =

mλ

? mλ ,

2

D?ν =

? ν mλ mλ

σ? σν

, (? = ν),

(4.13)

where D?ν are the correlation coe?cients between parameters λ? and λν . In order to accurately determine σλ a large number of simulations would be needed. If the measured values m λ obey Gaussian statistics then after Nt trials the variance is determined to a relative ac√ √ curacy 1/ N t and estimated values can di?er from their true values by σλ / N t . We have performed in excess of 5000 trials, for each input signal, and thus our results are accurate to better than 1 part in 70. Even more crucial than the number of simulations is the number of templates used and their range in the parameter space. We discuss these and other related issues next. The actual templates chosen, say for the parameter τ0 , in a given ‘experiment’ depend on the true parameters of the signal, the number of noise realisations employed and the expected value of the error. Let us suppose we have a ?rst guess of the error in τ0 , say στ0 . Then, we choose 51 uniformly spaced ?lters around τ0 (where τ0 is the signal chirp time) ? ? such that:

t τ0

? ∈ [?0 ? 5στ0 , τ0 + 5στ0 ] . τ

(4.14)

This implies that we are covering a 5σ width in τ0 at a resolution of στ /5. The probability, that a template between 4σ and 5σ from the true signal ‘clicks’, being ? 6 × 10?5 , we are on safe grounds since, in a given simulation, we consider no more than 5000 trials. (In comparison, the probability that a template between 3σ and 4σ clicks is 2.2 × 10?3 corresponding to an expected 13 events in 5000 trials.) For a post-Newtonian signal, which in e?ect needs to be spanned by a two-dimensional lattice of ?lters, the above choice of templates implies a requirement of 2601 × 2 ?lters in all. Here a factor of 2 arises because for each ?lter in the τ0 –τ1 space we will need two templates corresponding to the two independent values of the phase Φ: 0 and π/2. In the case of a Newtonian signal, the lattice being one-dimensional, one can a?ord a much higher resolution and range. Even with the aid of a mere 201 templates we can probe at a σ/10 resolution with a 10σ range. We start o? a simulation with the pretension that there is no knowledge of what the σλ ’s are. Thus, we choose as our ?rst trial a very large σλ and lay the lattice of templates. 18

With this lattice we perform a test run of 400 trials and examine the distribution of the measured values. If the distribution is not, as expected, a Gaussian then we alter σλ : we decrease it if the distribution is too narrow and increase it if the distribution is too wide and does not show the expected fallo?. In particular, we make sure that the templates at the boundary of the chosen range do not click even once and the skewness of the distribution is negligible. When for a certain σλ a rough Gaussian distribution is observed then we carry out a simulation with a larger number of trials (typically 5000). We subject the measured values in this larger simulation to the very same tests described above. We only consider for further analysis such simulations which ‘pass’ the above tests and determine the estimates, variances and covariances of the parameters using the measured values, with the aid of equation (4.13).

4. Numerical errors and remedies

There are several sources of numerical errors that tend to bias the results of a simulation unless proper care is exercised to rectify them. In this Section we point out the most important ones and show how they can be taken care of. Due to memory restrictions, the present version of our codes work with single precision except the FFT, which is implemented in double precision. In future implementations we plan to carry out all computations in double precision. This will possibly reduce some of the numerical noise that occurs, especially at high SNRs, in the present simulations. 1. Orthonormality of ?lters: For the sake of simplicity it is essential that the ?lters are normalised in the sense that their scalar product is equal to unity: q, q = 1. A waveform is normalized numerically using the discrete version of the scalar product: N = 1

N ?1 ?1 k=0 Sk

|?k | q

2

(4.15)

As mentioned earlier we use a two-dimensional basis of ?lters for the phase parameter. Choosing the two ?lters to be orthogonal to each other makes the maximisation over phase easier. However, here care must be exercised. Two ?lters q(t; ta , τ0 , τ1 , Φ = 0) and q(t; ta , τ0 , τ1 , Φ = π/2) are apparently orthogonal to each other. The numerically computed ‘angle’ between the two ?lters, chosen in this manner, often turns out to be greater than ? 10?2 radians. Consequently, one obtains erroneous correlations. In order to circumvent this problem we ?rst generate two ?lters that are roughly orthogonal to each other, as above, and then use the Gram-Schmidt method to orthogonalise the two vectors. If an explicit numerical orthogonalisation such as this is not implemented then the measured values of the various parameters show spurious oscillations in their distribution and the estimated values of the parameters tend to get biased. 2. Correlation function: The scalar product of two normalised templates q(t; t λ) and q(t; t λ′ ) is given by C(t λ? , t λ′ ) = q(t; t λ? ), q(t; t λ′ ) , ν ν q(t; t λ), q(t; t λ) = q(t; t λ′ ), q(t; t λ′ ) = 1, (4.16) where we have indicated the dependence of the scalar product on the various parameters by explicitly writing down the parameter subscripts. Let us ?x the parameters of one of the templates, say t λ, and vary the parameters of the other template. Of particular interest is the behaviour of C maximised over all but one of the parameters, say λν0 : Cmax (t λ? , t λ′ 0 ) = max C(t λ? , t λ′ ). ν ν ′

t λν

(4.17)

ν=ν0

19

Cmax is expected to drop monotonically as |λν0 ? λ′ 0 | increases. However, we have ν observed departures from such a behaviour possibly arising out of numerical noise. Such a behaviour causes bias in the estimation of parameters, and consequently in the determination of their covariances, especially at high SNRs. We have found no remedy to this problem and some of our results at high SNRs may have biases introduced by this e?ect. (Sampling the templates at a higher rate did not help in curing this problem.) 3. Grid e?ects: The parameters of a signal chosen for the purpose of simulation and detection can in principle be anything and in particular it need not correspond to any of the templates of the lattice. However, in practice we ?nd that whenever the signal parameters do not correspond to a member of the lattice then the resultant simulation has a bimodal distribution of the measured values. This is, of course, expected since a signal not on the grid is picked up by two nearest templates along each direction in the parameter space. Sometimes we do ?nd that the peaks corresponding to the bimodal distribution does not belong to the nearest neighbour ?lters but slightly away. This is related to the fact that the correlation function maximised over the time of arrival and the phase of the signal falls o? much too slowly along the τ0 –τ1 direction and small deviation from a monotonic fall can cause biases. (Such biases would be present in the case when a signal corresponds to one of the grid points though the magnitude of the e?ect would be lower.) In order to avoid this problem, and the consequent shifts in the estimation of parameters and errors in the determination of variances and covariances, we always choose the parameters of the signal to be that corresponding to a template. 4. Upper frequency cuto? and its e?ect on parameter estimation The Fisher information matrix computed using the stationary phase approximation in Section III does not include the e?ect of truncating the waveform at a = 6M — the plunge cuto?. As mentioned before, we have carried out simulations for both with, and without, incorporating the upper cuto?. As the covariance matrix incorporating the upper cuto? is not available we have been able to compare the Monte Carlo results with the covariance matrix only for the latter case, where the cuto? is held ?xed at 750 Hz. If we incorporate the upper cuto? into the Monte Carlo simulations the errors in the parameters are reduced drastically. The e?ect of the upper cuto? is expected to be more important for the higher mass binaries such as the ones we have considered. The ambiguity function, in this case, no longer remains independent of the point on the manifold. In other words, the correlation between two chirps depends not only on the di?erence between the parameters of the signals, but also on the absolute values of the parameters. The correlation surface also ceases to be symmetric i.e. the correlation between two chirps also depends on the sign of δλ, where, δλ is the di?erence in the values of the parameters. As the computational power required for carrying out simulations for lower mass binaries is not available to us the simulations have been restricted to ns-bh star binaries, where the e?ect of the upper cuto? is important. 5. Boundary e?ects For the purpose of simulations a grid of ?lters has to be set up ‘around’ the signal. The grid must be large enough so that the estimated parameters do not overshoot the boundary of the grid. This causes a problem as every value in the {τ0 , τ1 } plane does not lead to a meaningful value for the masses of the binary system. This does not however prevent us to construct a waveform with such a value for {τ0 , τ1 } even though the signal in general does not correspond to any ‘real’ binary system. This is valid, and even necessary, if we are to compare the numerical results with the covariance matrix. 6. Incorporating the cuto? in the presence of boundary e?ects If we wish to incorporate the e?ects of the upper cuto? in simulations then we run into a serious problem, as we would have to know the total mass of the binary in order to compute the upper cuto?. For an arbitrary {τ0 , τ1 } we can end up with negative and even complex values of the total mass and hence the upper cuto? at a = 6M is not 20

meaningful. Thus, we cannot even construct a waveform for an arbitrary combination of {τ0 , τ1 }. Therefore, in such cases, we restrict ourselves to simulations where the grid lies entirely within valid limits for {τ0 , τ1 }.

C. Results and Discussion

Our primary objective is to measure the variances and covariances following the method described in Sec. IV B 3 and study their departure from that predicted by analytical means (cf. Sec. IV A). We have carried out simulations for several values of the masses of the binary and in each case the signal strength (which is a measure of the SNR) is varied in the range 10–40. However, since the variances and covariances are independent of the absolute values of the parameters, for the parameter set that we employ, results are only quoted corresponding to a typical binary system. (See Sec. IV A for a discussion of the covariance matrix.) Similar results are obtained in other cases too. We use two sets of parameters to describe our results. Monte Carlo simulations allow us to directly measure the amplitude, the time of arrival, the phase at the time of arrival and the chirp time(s). This is the set {A, ta , Φ, τ0 , τ1 }. As we shall see below, the instant of coalescence can be measured much more accurately than the time of arrival. As a consequence of this, the direction to the source can be determined at a much greater accuracy by employing tC as a parameter instead of ta . Thus, we also quote estimates and errors for the parameter tC . Since the error in the estimation of the phase is quite large, even at high SNRs, we ignore it in our discussions. We ?rst deal with the Newtonian signal and highlight di?erent aspects of the simulation and discuss the results at length. We then consider the ?rst post-Newtonian corrected signal.

1. Newtonian signal

In the case of Newtonian signals the parameter space is e?ectively one-dimensional and, as mentioned earlier, in this case the lattice of templates covers a 10–σ range of the parameters at a resolution of σ/10 centred around the true parameters of the signal. In Fig. 5 we have shown the error σλ in the estimation of parameters ta , τ0 and tC , as a function of SNR, deduced using the covariance matrix as solid lines and computed using Monte Carlo simulations as dotted lines. The curve corresponding to the covariance matrix is obtained using an upper frequency cuto? fc = 750 Hz consistent with that used in our simulations. The errorbars in the estimation of σλ ’s are obtained using 4 simulations, each with 4000 trials. At low SNR’s σλ ’s have a larger uncertainty, as expected, and for ρ > 30 this uncertainty is negligible, and sometimes smaller than the thickness of the curves, except in the case of σtC (see below, for a possible explanation). At low SNRs (10–15) there is a large departure of the various σλ ’s from that inferred using the covariance matrix. At an SNR ? 17 the two curves merge (except in the case of σtC ) indicating the validity of the covariance matrix results for this and higher SNR’s. Interestingly, the agreement between Monte Carlo simulation results and those obtained using the covariance matrix is reached roughly at the same SNR irrespective of the parameter in question. We note that in spite of the fact that the time of arrival and the chirp time have large errors, the instant coalescence can be estimated very accurately—an order of magnitude better than either. What is puzzling, however, is that, in the case of tC , the Monte Carlo curve drops below the covariance matrix curve above an SNR of 15 and the two curves do not seem to converge to one another even at very high SNR’s. Coincident with the crossover of the two curves, the error in the estimation of σtC increases, contrary to what happens for the other parameters, signalling that there is a large ?uctuation in the estimation of σtC . This behaviour, we guess, is an artifact of the low value of the sampling rate. Of course, our sampling rate is su?ciently high to respect the sampling theorem. However, since tC is determined to an accuracy an order of magnitude better than either ta or τ0 , a much higher resolution in template-spacing would be needed for determining the error in the instant

21

of coalescence than that used for estimating the errors in the time of arrival or the chirp time(s). Testing this claim, unfortunately, is beyond the computer resources at our disposal since we would need a sampling rate of about 10 kHz with a ?lter-spacing 10?4 s. We hope to be able to resolve this issue in course of time. Nevertheless, the fact that the error in the estimation of σtC ?rst decreases with SNR and increases only after the two curves crossover, hints at the above possibility as a cause for this anomalous behaviour. This e?ect is also observed in the case of a post-Newtonian signal. ? In Table V we have given the actual signal parameters λ, estimated values of the parameters e λ (cf. equation (4.13)) and the corresponding errors in their estimation σλ , for several values of the SNR. Errors inferred from the covariance matrix can be read o? from Fig. 5. The estimated values are di?erent from the true values, some of them being overestimated and some others being underestimated. However, the deviations are often larger than what we expect. In a simulation that uses Nt trials the estimated parameters e λ, assuming a Gaussian distribution for the measured parameters m λ, can be di?erent from the true val√ ues by σλ / Nt . (In contrast, the measured values m λ can di?er from their true values by σλ or more.) However, we often obtain a slightly larger deviation, σλ < 2√ ? Nt

eλ

? ? σλ ? λ < 3√ , Nt

(4.18)

and we are unable to resolve this discrepancy. A more concrete test for the simulations is the histogram n(m λ) of the measured parameters, namely the frequency at which a given parameter clicks in a simulation. This is shown plotted in Fig. 6 for an SNR of 10. The skewness of the measured value is less than 10?2 . These results lend further support to the Monte Carlo simulations. There are visible asymmetries in the distributions of τ0 and ta and the asymmetries in the two cases are of opposite sense. This can, of course, be understood from the fact that ta and τ0 have a negative correlation coe?cient. The histogram of tC , even at an SNR of 10, has very few non-zero bins. This of course re?ects the fact that it is determined very accurately. We are unable to resolve the central peak in n(tC ) since, as mentioned earlier, the sampling rate and resolution in τ0 are not good enough to do so.

2. Post-Newtonian signal

As opposed to the Newtonian case here we have essentially a two-dimensional lattice of ?lters corresponding to τ0 and τ1 . For the purpose estimating variances and covariances we lay a mesh consisting of 2601 × 2 uniformly spaced ?lters around the true parameters of the signal. As pointed out in Sec. III B not all ?lters in the mesh, unlike in the Newtonian case, would correspond to the waveform from a realistic binary but that does not preclude their use in the Monte Carlo simulations. We shall see that the results of our simulations lend further support to the claim that for the purpose of detection, the parameter space need only be onedimensional [28]. The results obtained for the ?rst post-Newtonian signal are qualitatively similar to that of a Newtonian signal and we refer the reader, where appropriate, to the Newtonian case for a more complete discussion. In Fig. 7 we have shown the error in the estimation of parameters τ0 , τ1 , tC and ta , clockwise from top left, respectively, as a function of SNR. The solid and dotted curves are as in Fig. 5. Here again the upper frequency cuto? is taken to be 750 kHz. Just as in the case of a Newtonian signal here too the results obtained from Monte Carlo simulation are much higher than those obtained by employing the covariance matrix. At an SNR of 10, which is the expected value for a majority of events that initial LIGO and VIRGO interferometers will observe, the Monte Carlo values are more than thrice as much as their corresponding covariance matrix values and at an SNR of 15 the errors are roughly twice that expected from the covariance matrix. In absolute terms, however, the errors are still quite small compared to the actual parameter values: for a ns-ns binary, at an SNR of 10, στ0 στ1 ? 2.4%, ? 9.4%. (4.19) τ0 τ1 At an SNR of 10 the time of arrival can be measured to an accuracy of 72 ms in contrast to a value of 20 ms expected from the covariance matrix. As is well known, with the inclusion 22

of the post-Newtonian terms error in the estimation of the time of arrival and Newtonian chirp time increases by about a factor of 2 and 3, respectively [21,22]. As in the Newtonian case here again we see that the Monte Carlo curves approach the corresponding covariance matrix curves at a high SNR the only di?erence being that the agreement is reached at a higher SNR ? 25. For SNRs larger than this the two curves are in perfect agreement with each other. As mentioned earlier, σtC shows an anomalous behaviour possibly arising out of insu?cient resolution in the time of arrival and the chirp times. ? In Table V we have listed the true parameters λ, the estimated values e λ and the Monte Carlo errors σλ for di?erent SNRs. As in the Newtonian case here too the estimated values show a larger departure, than expected, from the true values. Histograms of the various measured parameters including tC are shown in Fig. 8 for a signal strength of 10. The skewness is below its standard deviation of 15/Nt [33] indicating the Gaussian nature of the various distributions. Even in the case of a post-Newtonian signal σtC is so small that we only have three non-zero bins in n(tC ). We now turn attention to other, more general, issues arising out of the simulations. In Sec. III C we have argued, on the basis of the behaviour of the noise-free correlation function, that the e?ective dimensionality of the parameter space for the purpose of detection, even in the case of a post-Newtonian signal, is only one-dimensional. The results of our Monte Carlo simulation unambiguously show that this is indeed true even in the presence of noise. We investigated the two-dimensional histogram, that gives the number of occurrances of di?erent templates in the lattice, in a particular simulation. The templates that ‘click’ are all aligned along the line τ0 + τ1 = const. In a total of 5000 realisations corresponding to this simulation there is only one instance when a ?lter outside this region clicks giving a probability of less than 10?3 for a template outside this region to give a maximum. Consequently, it is only necessary to choose a single ?lter along this curve. ? The distribution of the maximum correlation Cmax (λ, t λ) obtained from di?erent noise realisations needs a special mention since it has an inherent bias. In Fig. 9 we have shown the distribution of the maximum correlation taken from one of our simulations corresponding to an SNR of 10. Notice a slight shift of the distribution towards a higher value and this cannot be accommodated within the expected ?uctuation in the mean. The measured value of the standard deviation σA is 0.95. Since the number of simulations is 5000 we expect that the √ signal strength should di?er from the true value of 10 by no more than σA / 5000 = 0.014. However, the mean value is 10.26 giving a deviation of 0.26 which is about 20 times larger than that expected. This occurs at all SNRs and for both Newtonian and post-Newtonian signals. This of course does not mean there is bug in the way we are computing the maximum correlation. In the process of maximisation, values greater than the signal strength are favoured and consequently the mean of the maximum correlation shows a shift towards a higher value. This suggests that that the maximum of the correlation is a biased estimator of the signal strength. We ?nd, consistent with the covariance matrix calculation, that the amplitude parameter is uncorrelated with the rest of the parameters; crosscorrealtion coe?cients D0 ?, ? = 0, (cf. eq. 4.2) inferred from our Monte Carlo simulations are less than ? 10?6 . Finally, it is of interest to note how the phase parameter Φ is correlated with the time of arrival. A plot of m Φ versus m ta is shown in Fig. 10. We ?nd that the measured values of the time of arrival and the phase are such that 2πf0 m ta = m Φ, where f0 has a value of approximately 51 Hz. When the time of arrival shifts by more than a cycle of the signal the phase jumps by a factor 2π leading to the points seen in the top-left and bottom-right corner of the ?gure. This makes the estimation of the phase and the error in its estimation pretty involved.

3. Incorporating the e?ects of upper cuto?

As mentioned before, incorporating an upper cuto? at the onset of the plunge has a drastic e?ect on the estimation of parameters. The incorporation of the upper cuto? is implemented by stopping the waveform when the instantaneous frequency reaches the fre-

23

quency associated with the onset of the plunge or 750 Hz, whichever is lower. However, due to computational constraints we have carried out the simulations only for high mass binaries and hence the upper cuto? plays an important role in all our simulations. It is to be noted that the discussion of the ambiguity function in section III C is not valid when the e?ect of the upper cuto? is taken into account though for low mass binaries, such as ns-ns binaries, the results there will still be valid. The further dependence of the signal waveform on the total mass of the system through the upper cuto? means that we can estimate the individual masses more exactly though the computational power is bound to increase. In order to carry out the simulations for the present case we selected a 10M⊙ ? 1.2M⊙ binary system as this enables us to choose the ?lter grid well within valid limits of τ0 and τ1 . The simulations were carried out for various values of SNR starting from 10. The histograms of the estimated parameters at an SNR of 10 is shown in ?g. 11. At this SNR the errors obtained are στ0 = 39.3 ms, στ1 = 22.4 ms, σta = 23.1 ms and σtC = 0.6 ms. These can be compared with the values in table V and we can see that except for the parameter tC the errors are substantially lesser when the upper cuto? is incorporated into the waveform. It is necessary to recompute the covariance matrix, as emphasized before, including the e?ect of the upper cuto? in order to compare these numerically obtained values with the covariance matrix. In order to do this it is not enough to replace the upper limit in the integral in eq. (3.19) with the upper cuto?. The waveform now depends on the total mass of the system through the upper cuto? and this information has to be incorporated into the waveform. We carried out the simulations for various SNRs for the same value of masses quoted above. In the absence of the estimates of the covariance matrix when the upper cuto? is incorporated we assume that at an SNR of 40 the Monte Carlo estimates are consistent with those of the covariance matrix. In Figure 12 we illustrate the results of our simulations. The points on the dotted line are the values obtained through Monte Carlo estimates and the continuous line is obtained by ?tting a 1/ρ dependence of the errors on the SNR assuming consistency at an SNR of 40. It is seen as in the previous simulations that except for the parameter tC the other parameters are fairly consistent with a 1/ρ dependence of the errors on the SNR at an SNR greater than 25.

V. CONCLUSIONS

In this paper we have introduced the use of di?erential geometry in studying signal analysis and have addressed issues pertaining to optimal detection strategies of the chirp waveform. We have also carried out Monte Carlo simulations to check how well the covariance matrix estimates the errors in the parameters of the chirp waveform. We summarize below our main results. 1. We have developed the concept of a signal manifold as a subset of a ?nite dimensional vector space of detector outputs. Using the correlation between two signal vectors as a scalar product we have induced a metric on the signal manifold. With this geometric picture it is possible to pose the question of optimal detection in a more general setting. We suggest that the set of template waveforms for the detection of the chirp signal need not correspond to any point on the chirp manifold. We propose an algorithm to choose templates o? the signal manifold and show that the drop in the correlation due to the discreteness of the set of templates is reduced. This algorithm, though certainly not the best, motivates the search of more e?cient templates. In addition, the chirp manifold corresponding to the second post-Newtonian waveform is shown to be e?ectively one-dimensional. This has important implications for the computational requirement for the on-line detection of the chirp signal. The use of a convenient set of parameters of the chirp waveform for carrying out numerical and analytical simulations is stressed. These parameters are such that the metric components are independent of the parameters which implies that the manifold is ?at and the corresponding ‘coordinate system’ is Cartesian. As the metric de?ned is nothing but the Fisher information matrix, the covariance matrix, being the inverse of the Fisher information matrix, is also independent of the parameters.

24

2. Monte Carlo simulations have been carried out for the case of the initial LIGO to ?nd out whether the actual errors in the estimation of parameters is consistent with the values predicted by the covariance matrix. Simulations have been carried out for both the Newtonian as well as the post-Newtonian waveforms. We have restricted ourselves to the case of high mass binary systems, such as bh-ns binaries, where the computational requirement is not very heavy since the length of the data train, in such cases, works out to be less than 8 s. Nevertheless, as has been shown in this paper, the covariance matrix is independent of the parameters identi?ed by us when waveforms are terminated at a constant upper cuto? irrespective of their masses. Consequently, our results will hold good for binary systems of arbitrary masses. We point out the major problems that arise while performing a numerical simulation and, where appropriate, we suggest how they may be taken care of. In particular, the e?ect of incorporating the upper cuto? in the frequency of the gravitational wave at the onset of the plunge, which essentially depends on the total mass of the binary, is extremely important for high mass binaries. Since the covariance matrix with the inclusion of such a mass-dependent upper cuto? is not available we have carried out most of our simulations using a constant upper cuto?. This enables us to directly compare the results of our Monte Carlo simulations with those of the analytically computed covariance matrix. Since for binaries with total mass less than 5M⊙ the plunge induced upper cuto? is larger than that induced by the detector noise these e?ects can be ignored for such binaries. The numerical experiments indicate that the covariance matrix underestimates the errors in the determination of the parameters even at SNRs as high as 20. In the Newtonian case the correlation coe?cient of the time of arrival ta and the Newtonian chirp time τ0 is found to be very close to ?1, so much so that even at a SNR of 7.5, the instant of coalescence tC = ta + τ0 remains practically a constant. The error in the estimation of τ0 for the post-Newtonian waveform is about four times the error obtained in the case of the Newtonian waveform at the same SNR. This is expected as the ?rst post-Newtonian correction to the waveform introduces a new parameter τ1 (called the ?rst post-Newtonian chirp time) which is highly (anti)correlated with τ0 . For the post-Newtonian waveform at an SNR of 10 the error in τ0 is about 3 times that predicted by the covariance matrix. This corresponds to a factor of 2 in the chirp mass M. The distributions for the parameters have been obtained and are seen to be unimodal distributions and are slightly more sharper than a Gaussian. When the plunge induced upper cuto? is incorporated into the waveform the errors in the estimation of parameters decrease by a factor of about 2.5. The correlation coe?cient between τ0 and τ1 is also found to decrease, which is consistent with our discussion in Section IV A. The results obtained suggest that higher moments than that used in obtained in computing the covariance matrix may be important in the determination of the errors in parameter estimation. In the geometric picture this amounts to taking into account curvature e?ects, either intrinsic or extrinsic. 3. We suggest that tC is a more suitable parameter to estimate the direction to the source than the time of arrival. The latter is a kinematical parameter that ?xes the time at which the gravitational wave frequency reaches the lower cuto? of the detector while the parameter tC has the physical signi?cance of being the instant of coalescence. At an SNR of 10 the error in ta is too large (20 ms) to deduce the direction to the source accurately, whereas the error in the parameter tC is less than 0.5 ms. This will further go down substantially for the advanced LIGO. A detailed analysis of the determination of the direction to the binary using delays in tC between di?erent detectors is carried out in Bhawal and Dhurandhar [35] (also see [29]). We now suggest further work which needs to be done along the lines of this paper. A full understanding of the chirp signal manifold when higher post-Newtonian corrections are incorporated into the waveform is in order. This will help in the development of more e?cient algorithms for the choice of templates in the detection problem and facilitate reduction in

25

computational time. The Monte Carlo simulations which we have carried our are for the case of a binary waveform correct up to ?rst post-Newtonian order. Moreover, only circular orbits are considered. The e?ect of eccentricity is currently being investigated [36]. Performing simulations when higher post-Newtonian corrections are taken into account calls for an immense amount of computational time. Fortunately, matched ?ltering algorithm being amenable to parallelization [37], one could aim at using the massively parallel computers, that are now becoming available the world over, in performing such simulations.

ACKNOWLEDGMENTS

The authors would like to thank the members of the gravitational wave group at IUCAA especially S.D. Mohanty and B. Bhawal for many helpful discussions. R.B. is being supported by the Senior Research Fellowship of CSIR, India.

[1] A. Abramovici et. al., Science 256 325 (1992). [2] C. Bradaschia et. al., Nucl. Instrum. Methods Phys. Res., Sect A 518 (1990). [3] B.F. Schutz, in Gravitational Collapse and Relativity, Edited by H. Sato and T. Nakamura, (World Scienti?c, Singapore, 1986), pp. 350-368. [4] D. Markovi?, Phys. Rev. D 48, 4738 (1993). c [5] L. Blanchet and B.S. Sathyaprakash, Phys. Rev. Lett., 74, 1067 (1995). [6] K.S. Thorne, Gravitational waves from compact objects, to be published in Proceedings of IAU symposium 165, Compact stars in binaries, edited by J. van Paradijs, E. van den Heuvel, and E. Kuulkers, (Kluwer Academic Publishers). [7] K.S. Thorne, Gravitational waves, to be published in Proceedings of the Snowmass 95 Summer Study on Particle and Nuclear Astrophysics and Cosmology, edited by E.W. Kolb and R. Peccei, (World Scienti?c, Singapore). [8] S.D. Mohanty and B.S Sathyaprakash, A modi?ed periodogram for the detection of gravitational waves from coalescing binaries, submitted to Phys. Rev. D (1995). [9] S, Smith, Phys. Rev. D 36, 2901 (1987). [10] K.S. Thorne, in 300 Years of Gravitation, S.W. Hawking and W.Israel (eds.), (Cambridge Univ. Press, 1987). [11] C.W. Helstrom, Statistical Theory of Signal Detection, 2nd. ed, (Pergamon Press, London, 1968). [12] B.F. Schutz, in The Detection of Gravitational Radiation, edited by D. Blair (Cambridge, 1989) pp 406-427. [13] B.S. Sathyaprakash and S.V. Dhurandhar, Phys. Rev. D 44, 3819 (1991). [14] S.V. Dhurandhar and B.S. Sathyaprakash, Phys. Rev. D 49, 1707 (1994). [15] R. Balasubramanian and S.V. Dhurandhar, Phys. Rev. D 50, 6080 (1994). [16] K. Kokkotas, A. Krolak, and G. Tsegas, Class. Quant. Grav., 11, 1901 (1994). [17] T. A. Apostolatos, Phys. Rev. D, 52, 605 (1995). [18] A. Krolak, K.D. Kokkotas and G. Sch¨fer, Phys. Rev. D. in press. a [19] Shun-ichi Amari, Di?erential Geometric Methods in Statistics, (Springer-Verlag). [20] L. S. Finn, Phys. Rev. D 46, 5236 (1992). [21] L.S. Finn and D.F. Cherno?, Phys. Rev. D 47, 2198 (1993). [22] C. Cutler and E. Flanagan, Phys.Rev. D 49, 2658 (1994). [23] L. Blanchet and B.S. Sathyaprakash, Class. Quantum Grav., 11, 2807 (1994). [24] A. Krolak, in Gravitational wave data analysis, edited by B.F. Schutz, (Dordrecht : Kluwer, 1989), pp. 59-69. [25] A. Krolak, J. A. Lobo and B. J. Meers, Phys. Rev D 48, 3451 (1993). [26] E. Poisson and C.M. Will, Phys. Rev. D 52, 848 (1995). [27] P. Jaranowski, K.D. Kokkotas, A. Krolak and G. Tsegas, report (unpublished). [28] B.S. Sathyaprakash, Phys. Rev. D 50, R7111 (1994).

26

[29] R. Balasubramanian, B.S. Sathyaprakash and S.V. Dhurandhar, submitted for publication to, Pramana J. Phys. (1995). [30] C. Cutler et al, Phys. Rev. Lett. 70, 2984 (1993). [31] L. Blanchet, T. Damour and B.R. Iyer, Phys. Rev. D 51, 5360 (1995). [32] S. V. Dhurandhar and B.F. Schutz, Phys. Rev. D 50, 2390 (1994) [33] William H. Press, Saul A. Teukolsky, William T. Vellerling and Brian P. Flannery Numerical Recipes (Cambridge University Press, 1993). [34] L. Blanchet, T. Damour, B.R. Iyer, C.M. Will and A.G. Wiseman, Phys. Rev. Lett. 74, 3515, (1995). [35] B. Bhawal and S.V. Dhurandhar, private communication. [36] R. Balasubramanian, A. Gopakumar, B.R. Iyer and B.S. Sathyaprakash, in preparation. [37] B.S. Sathyaprakash and S.V. Dhurandhar, J. Computational Phys., 109, 215 (1993).

TABLE I. The estimated value of the parameters and their errors e λ (σλ ) for the Newtonian ? waveform. The actual values of the parameters taken were τ0 = 5558.0 ms and ta = 200.0 ms. ? Except for the parameter A all values are quoted in ms. ρ = 7.5 (στ0 ) (σta ) e tC (σtC ) e A (σA )

e τ0 e ta

ρ = 10.0 5557.3 (15.6) 200.7 (15.5) 5758.0 (0.20) 10.12 (0.98)

ρ = 12.5 5557.6 (9.0) 200.4 (8.9) 5758.0 (0.10) 12.582 (0.99)

ρ = 15.0 5557.7 (6.2) 200.3 (6.2) 5758.0 (0.08) 15.067 (0.99)

ρ = 20.0 5558.5 (3.9) 199.5 (3.9) 5758 (0.02) 20.05 (0.99)

5557.4 (29.7) 200.6 (29.4) 5758.0 (0.30) 7.696 (0.96)

TABLE II. The estimated value of the parameters and the errors in their estimation e λ (σλ ) for the post-Newtonian waveform. The actual values of the parameters taken were τ0 = 5558.0 ms, τ1 = 684.0 ms and ta = 300.0 ms Except for the parameter A all values are quoted in ms. ρ = 10.0 (στ0 ) (στ1 ) e ta (σta ) e tC (σtC ) e A (σA )

e τ0 e τ1

ρ = 15.0 5555.9 (64.8) 685.2 (32.6) 300.88 (33.4) 6.542 (0.30) 15.15 (0.98)

ρ = 20.0 5557.2 (30.1) 684.6 (16.4) 300.28 (14.5) 6.542 (0.17) 20.11 (.99)

ρ = 25.0 5556.7 (19.3) 684.8 (10.5) 300.5 (9.4) 6.542 (0.10) 20.1 (0.99)

ρ = 30.0 5558.7 (13.6) 683.6 (7.2) 299.7 (6.8) 6.542 (0.06) 30.1 (0.99)

5554.9 (136.1) 685.55 (65) 301.49 (73.26) 6.542 (0.54) 10.25 (0.97)

27

FIG. 1. This ?gure illustrates the correlation function along the line of curvature as a function of τ0 , the Newtonian chirp time, for the following three cases: (i) when the ?lter is placed on the manifold (dotted line), (ii) when the ?lter is the ‘average’ signal vector over the region (dashed line) and (iii) when the ?lter is chosen to be an appropriate linear combination of the previous two ?lters (solid line).

FIG. 2. Contour diagram of the ambiguity function for the second post-Newtonian case.

FIG. 3. Dependence of the errors in the estimation of the parameters of the post-Newtonian waveform on the upper cuto? frequency. The SNR is kept ?xed at a value of 10.

FIG. 4. Dependence of the errors in the estimation of the parameters of the post-Newtonian waveform on the total mass of the binary keeping the distance to the binary ?xed. The waveforms are cuto? at frequencies corresponding to the onset of the plunge orbit and the SNR is normalized at a value of 10 for a 1.4 M⊙ -10 M⊙ binary.

FIG. 5. Dependence of the errors in the estimation of parameters of the Newtonian waveform i.e. {στ0 , σta , σtC } as a function of SNR. The solid line represents the analytically computed errors whereas the dotted line represents the errors obtained through Monte Carlo simulations.

FIG. 6. Distributions of the measured values of the parameters for the case of Newtonian signal. The total number of noise realizations is 5000. FIG. 7. Dependence of the errors in the estimation of parameters of the post-Newtonian waveform i.e. {στ0 , στ1 , σta , σtC } as a function of SNR. The solid line represents the analytically computed errors whereas the dotted line represents the errors obtained through Monte Carlo simulations.

FIG. 8. Distributions of the measured parameters of the post-Newtonian waveform at a SNR of 25. The number of noise realizations is 5000. FIG. 9. The histogram of random variable m A at an SNR of 20. The variance in this parameter is independent of the SNR and is approximately equal to unity.

FIG. 10. The correlation between ta and Φ is illustrated. The phase parameter simply follows the time of arrival parameter.

FIG. 11. Distributions for the measured parameters of the post-Newtonian waveform incorporating the e?ect of the upper cuto?. The errors in the determination of the parameters is much smaller in this case. The total number of noise realizations is 5000.

FIG. 12. The dependence of the errors on the SNR (ρ) when the upper cuto? is incorporated into the waveform

28

1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

赞助商链接

- Wave Effects in Gravitational Lensing of Gravitational Waves from Chirping Binaries
- Hierarchical search strategy for the detection of gravitational waves from coalescing binar
- Localizing coalescing massive black hole binaries with gravitational waves
- Detecting gravitational waves from precessing binaries of spinning compact objects Adiabati
- Phasing of gravitational waves from inspiralling eccentric binaries at the third-and-a-half

更多相关文章：
更多相关标签：