Gravitational-Wave Data Analysis. Formalism and Sample Applications: The Gaussian Case

Abstract

The article reviews the statistical theory of signal detection in application to analysis of deterministic gravitational-wave signals in the noise of a detector. Statistical foundations for the theory of signal detection and parameter estimation are presented. Several tools needed for both theoretical evaluation of the optimal data analysis methods and for their practical implementation are introduced. They include optimal signal-to-noise ratio, Fisher matrix, false alarm and detection probabilities, ${\mathcal F}$-statistic, template placement, and fitting factor. These tools apply to the case of signals buried in a stationary and Gaussian noise. Algorithms to efficiently implement the optimal data analysis techniques are discussed. Formulas are given for a general gravitational-wave signal that includes as special cases most of the deterministic signals of interest.

1 Introduction

In this review we consider the problem of detection of deterministic gravitational-wave signals in the noise of a detector and the question of estimation of their parameters. The examples of deterministic signals are gravitational waves from rotating neutron stars, coalescing compact binaries, and supernova explosions. The case of detection of stochastic gravitational-wave signals in the noise of a detector is reviewed in [8]. A very powerful method to detect a signal in noise that is optimal by several criteria consists of correlating the data with the template that is matched to the expected signal. This matched-filtering technique is a special case of the maximum likelihood detection method. In this review we describe the theoretical foundation of the method and we show how it can be applied to the case of a very general deterministic gravitational-wave signal buried in a stationary and Gaussian noise.

Early gravitational-wave data analysis was concerned with the detection of bursts originating from supernova explosions [144]. It involved analysis of the coincidences among the detectors [70]. With the growing interest in laser interferometric gravitational-wave detectors that are broadband it was realized that sources other than supernovae can also be detectable [135] and that they can provide a wealth of astrophysical information [122, 77]. For example, the analytic form of the gravitational-wave signal produced during the inspiral phase of a compact binary coalescence is known in terms of a few parameters to a good approximation (see, e.g., [30] and Section 2.4 of [66]). Consequently one can detect such a signal by correlating the data with the predicted waveform (often called the template) and maximizing the correlation with respect to the parameters of the waveform. Using this method one can pick up a weak signal from the noise by building a large signal-to-noise ratio over a wide bandwidth of the detector [135]. This observation has led to rapid development of the theory of gravitational-wave data analysis. It became clear that the detectability of sources is determined by optimal signal-to-noise ratio, which is the power spectrum of the signal divided by the power spectrum of the noise integrated over the bandwidth of the detector.

An important landmark was a workshop entitled Gravitational Wave Data Analysis held in Dyffryn House and Gardens, St. Nicholas near Cardiff, in July 1987 [123]. The meeting acquainted physicists interested in analyzing gravitational-wave data with the basics of the statistical theory of signal detection and its application to detection of gravitational-wave sources. As a result of subsequent studies, the Fisher information matrix was introduced to the theory of the analysis of gravitational-wave data [50, 76]. The diagonal elements of the Fisher matrix give lower bounds on the variances of the estimators of the parameters of the signal and can be used to assess the quality of astrophysical information that can be obtained from detections of gravitational-wave signals [41, 75, 25]. It was also realized that the application of matched-filtering to some sources, notably to continuous sources originating from neutron stars, will require extraordinary large computing resources. This gave a further stimulus to the development of optimal and efficient algorithms and data analysis methods [124].

A very important development was the work by Cutler et al. [43] where it was realized that for the case of coalescing binaries matched filtering was sensitive to very small post-Newtonian effects of the waveform. Thus, these effects can be detected. This leads to a much better verification of Einstein’s theory of relativity and provides a wealth of astrophysical information that would make a laser interferometric gravitational-wave detector a true astronomical observatory complementary to those utilizing the electromagnetic spectrum. As further development of the theory, methods were introduced to calculate the quality of suboptimal filters [13], to calculate the number of templates required to do a search using matched-filtering [102], to determine the accuracy of templates required [33], and to calculate the false alarm probability and thresholds [68]. An important point is the reduction of the number of parameters that one needs to search for in order to detect a signal. Namely estimators of a certain type of parameters, called extrinsic parameters, can be found in a closed analytic form and consequently eliminated from the search. Thus, a computationally-intensive search need only be performed over a reduced set of intrinsic parameters [76, 68, 78].

Techniques reviewed in this paper have been used in the data analysis of prototypes of gravitational-wave detectors [100, 99, 12] and in the data analysis of gravitational-wave detectors currently in operation [133, 23, 4, 3, 2].

2 Response of a Detector to a Gravitational Wave

There are two main methods of detecting gravitational waves currently in use. One method is to measure changes induced by gravitational waves on the distances between freely-moving test masses using coherent trains of electromagnetic waves. The other method is to measure the deformation of large masses at their resonance frequencies induced by gravitational waves. The first idea is realized in laser interferometric detectors (both Earth-based [105, 141, 53] and space-borne [85, 137] antennas) and Doppler tracking experiments [14], whereas the second idea is implemented in resonant mass detectors (see, e.g., [22]).

2.1 Doppler shift between freely falling particles

We start by describing the change of a photon’s frequency caused by a passing gravitational wave and registered by particles (representing different parts of a gravitational-wave detector) freely falling in the field of the gravitational wave. The detailed derivation of the formulae we show here can be found in Chapter 5 of [66] (see also [47, 15, 120]). An equivalent derivation of the response of test masses to gravitational waves in the local Lorentz gauge (without making use of the long-wavelength approximation) is given in [113].

We employ here the transverse traceless (TT) coordinate system (more about the TT gauge can be found, e.g., in Section 35.4 of [91] or in Section 1.3 of [66]). A spacetime metric describing a plane gravitational wave traveling in the +z direction of the TT coordinate system (with coordinates x⁰ ≡ ct, x¹ ≡ x, x² ≡ y, x³ ≡ z), is described by the line element

$${\rm{d}}{{\rm{s}}^2}= - {c^2}{\rm{d}}{t^2} + \left({1 + {h_ +}\left({t - {z \over c}} \right)} \right){\rm{d}}{x^2} + \left({1 - {h_ +}\left({t - {z \over c}} \right)} \right){\rm{d}}{y^2} + 2{h_ \times}\left({t - {z \over c}} \right){\rm{d}}x{\rm{d}}y + {\rm{d}}{z^2},$$

(1)

where h₊ and h_× are the two independent polarizations of the wave. We assume that the wave is weak, i.e., for any instant of time t,

$$\vert {h_ +}(t)\vert \ll 1,\,\vert {h_ \times}(t)\vert \ll 1.$$

(2)

We will neglect all terms of order h² or higher. The form of the line element (1) implies that the functions h₊(t) and h_×(t) describe the wave-induced perturbation of the flat Minkowski metric at the origin of the TT coordinate system (where x = y = z = 0). It is convenient to introduce the three-dimensional matrix of the spatial metric perturbation produced by the gravitational wave (at the coordinate system’s origin),

$${\rm{H}}(t): = \left({\begin{array}{*{20}c} {{h_ +}(t)} & {{h_ \times}(t)} & 0 \\ {{h_ \times}(t)} & {- {h_ +}(t)} & 0 \\ 0 & 0 & 0 \end{array}} \right).$$

(3)

Let two particles freely fall in the field (1) of the gravitational wave, and let their spatial coordinates remain constant, so the particles’ world lines are described by equations

$$t({\tau _a}) = {\tau _a},\,x({\tau _a}) = {x_a},\,y({\tau _a}) = {y_a},\,z({\tau _a}) = {z_a},\,a = 1,2,$$

(4)

where (x_a, y_a, z_a) are spatial coordinates of the ath particle and τ_a is its proper time. These two particles measure, in their proper reference frames, the frequency of the same photon traveling along a null geodesic x^α = x^α(λ), where λ is some affine parameter. The coordinate time, at which the photon’s frequency is measured by the ath particle, is equal to t_a (a = 1, 2); we assume that t₂ > t₁. Let us introduce the coordinate time duration t₁₂ of the photon’s trip and the Euclidean coordinate distance L₁₂ between the particles:

$${t_{12}}: = {t_2} - {t_1},{L_{12}}: = \sqrt {{{({x_2} - {x_1})}^2} + {{({y_2} - {y_1})}^2} + {{({z_2} - {z_1})}^2}.}$$

(5)

Let us also introduce the 3-vector n of unit Euclidean length directed along the line connecting the two particles. We arrange the components of this vector into the column 3 × 1 matrix n (thus, we distinguish here the 3-vector n from its components being the elements of the matrix n; the same 3-vector can be decomposed into components in different spatial coordinate systems):

$${\rm{n: = (cos}}\,\alpha {\rm{,cos}}\beta {\rm{,cos}}\gamma)^{\rm{T}} = \left({\begin{array}{*{20}c} {\cos \alpha} \\ {\cos \beta} \\ {\cos \gamma} \\ \end{array}} \right),$$

(6)

where the superscript T denotes matrix transposition. If one neglects the spacetime curvature caused by the gravitational wave, then α, β, and γ are the angles between the path of the photon in the 3-space and the coordinate axis x, y, or z, respectively (obviously, α, β, γ ∈ 〈0;; π〉 and cos² α + cos² β + cos² γ = 1). Let us denote the value of the frequency registered by the ath particle by ν_a (a = 1, 2) and let us finally define the relative change of the photon’s frequencies,

$${y_{12}}: = {{{\nu _2}} \over {{\nu _1}}} - 1.$$

(7)

then, it can be shown (see Chapter 5 of [66] for details) that the frequency ratio y₁₂ can be written [making use of the quantities introduced in Eqs. (3) and (5)–(6)] as follows (the dot means here matrix multiplication):

$${y_{12}} = {1 \over {2(1 - \cos \gamma)}}{{\rm{n}}^{\rm{T}}} \cdot \left({\rm{H}\left({{t_1} - {{{z_1}} \over c}} \right) - {\rm{H}}\left({{t_1} - {{{z_1}} \over c} + (1 - \cos \gamma){{{L_{12}}} \over c}} \right)} \right) \cdot {\rm{n + \mathcal{O}(}}{{\rm{h}}^2}).$$

(8)

It is convenient to introduce the unit 3-vector k directed from the origin of the coordinate system to the source of the gravitational wave. In the coordinate system adopted by us the wave is traveling in the +z direction. Therefore, the components of the 3-vector k, arranged into the column matrix k, are

$${\rm{k}} = {(0,0, - 1)^{\rm{T}}}.$$

(9)

The positions of the particles with respect to the origin of the coordinate system we describe by the 3-vectors x_a (a = 1, 2), the components of which we put into the column matrices x_a:

$${{\rm x} _a} = {({x_a},{y_a},{z_a})^{\rm{T}}},\,a = 1,2.$$

(10)

Making use of Eqs. (9)–(10) we rewrite the basic formula (8) in the following form

$${y_{12}} = {{{{\rm{n}}^{\rm{T}}} \cdot \left({{\rm{H}}\left({{t_1} + {{{{\rm{k}}^{\rm{T}}} \cdot {{\rm{x}}_1}} \over c}} \right) - {\rm{H}}\left({{t_1} + {{{L_{12}}} \over c} + {{{{\rm{k}}^{\rm{T}}} \cdot {{\rm{x}}_2}} \over c}} \right)} \right) \cdot {\rm{n}}} \over {2(1 + {{\rm{k}}^{\rm{T}}} \cdot {\rm{n}})}} + {\mathcal O}({h^2}).$$

(11)

To obtain the response for all currently working and planned detectors it is enough to consider a configuration of three particles shown in Figure 1. Two particles model a Doppler tracking experiment, where one particle is the Earth and the other is a distant spacecraft. Three particles model a ground-based laser interferometer, where the masses are suspended from seismically-isolated supports or a space-borne interferometer, where the three test masses are shielded in satellites driven by drag-free control systems. In Figure 1 we have introduced the following notation: O denotes the origin of the TT coordinate system related to the passing gravitational wave, x_a(a = 1, 2, 3) are 3-vectors joining O and the particles, n_a and L_a (a = 1, 2, 3) are, respectively, 3-vectors of unit Euclidean length along the lines joining the particles and the coordinate Euclidean distances between the particles, where a is the label of the opposite particle. We still assume that the spatial coordinates of the particles do not change in time.

Figure 1

Schematic configuration of three freely-falling particles as a detector of gravitational waves. The particles are labelled 1, 2, and 3, their positions with respect to the origin O of the coordinate system are given by 3-vectors x_a (a = 1, 2, 3). The Euclidean separations between the particles are denoted by L_a, where the index a corresponds to the opposite particle. The unit 3-vectors n_a point between pairs of particles, with the orientation indicated.

Let us denote by ν₀ the frequency of the coherent beam used in the detector (laser light in the case of an interferometer and radio waves in the case of Doppler tracking). Let the particle 1 emit the photon with frequency ν₀ at the moment t₀ towards the particle 2, which registers the photon with frequency ν′ at the moment ${t\prime} = {t_0} + {L_3}/c + {\mathcal O(h)}$. The photon is immediately transponded (without change of frequency) back to the particle 1, which registers the photon with frequency ν at the moment $t = {t_0} + 2{L_3}/c + {\mathcal O(h)}$. We express the relative changes of the photon’s frequency y₁₂ ≔ (ν′ − ν₀)/ν₀ and y₂₁ ≔ (ν − ν′)/ν′ as functions of the instant of time t. Making use of Eq. (11) we obtain

$${y_{12}}(t) = {1 \over {2(1 - {{\rm{k}}^{\rm{T}}} \cdot {{\rm{n}}_3})}}{\rm{n}}_3^{\rm{T}} \cdot \left({{\rm{H}}\left({t - {{2{L_3}} \over c} + {{{{\rm{k}}^{\rm{T}}} \cdot {{\rm{x}}_1}} \over c}} \right) - {\rm{H}}\left({t - {{{L_3}} \over c} + {{{{\rm{k}}^{\rm{T}}} \cdot {{\rm{x}}_2}} \over c}} \right)} \right) \cdot {{\rm{n}}_3} + \mathcal{O}({h^2}),$$

(12A)

$${y_{21}}(t) = {1 \over {2(1 + {{\rm{k}}^{\rm{T}}}\cdot{{\rm{n}}_3})}}{\rm{n}}_3^{\rm{T}}\cdot\left({{\rm{H}}\left({t - {{{L_3}} \over c} + {{{{\rm{k}}^{\rm{T}}}\cdot{{\rm{x}}_2}} \over c}} \right) - {\rm{H}}\left({t + {{{{\rm{k}}^{\rm{T}}}\cdot{{\rm{x}}_1}} \over c}} \right)} \right)\cdot{{\rm{n}}_3} + {\mathcal O}({h^2}).$$

(12B)

The total frequency shift y₁₂₁ ≔ (ν − ν₀)/ν₀ of the photon during its round trip can be computed from the one-way frequency shifts y₁₂ and y₂₁ given above:

$${y_{121}} = {\nu \over {{\nu _0}}} - 1 = {\nu \over {{\nu {\prime}}}}{{{\nu {\prime}}} \over {{\nu _0}}} - 1 = ({y_{21}} + 1)({y_{12}} + 1) - 1 = {y_{12}} + {y_{21}} + {\mathcal O}({h^2}).$$

(13)

2.2 Long-wavelength approximation

Let L be the size of the detector and ƛ ≔ λ/(2π) be the reduced wavelength of the gravitational wave impinging on the detector. In the long-wavelength approximation the condition ƛ ≫ L is fulfilled. The angular frequency of the wave equals ω = c/ƛ. Time delays caused by the finite speed of the wave propagating across the detector are of order Δt ∼ L/c, but

$$\omega \Delta t \sim {L \over -{{\rlap{-}\lambda}}} \ll 1,$$

(14)

so time delays across the detector are much shorter than the period of the gravitational wave and can be neglected. It means that with a good accuracy the gravitational-wave field can be treated as being uniform (but time-dependent) in the space region that covers the entire detector. To detect gravitational waves with some dominant angular frequency ω one must collect data over time intervals longer (sometimes much longer) than the gravitational-wave period. This implies that in Eq. (8) for the relative frequency shift, the typical value of the quantity $\bar t: = {t_1} - {z_1}/c$ will be much larger than the retardation time Δt ≔ L₁₂/c Therefore, we can expand this equation with respect to Δt and keep terms only linear in Δt. After doing this one obtains (see Section 5.3 in [66] for more details):

$${y_{12}} = - {{{L_{12}}} \over {2c}}{{\rm{n}}^{\rm{T}}}\cdot{\rm{\dot H}}({t_1} - {{{z_1}} \over c})\cdot{\rm{n}} + {\mathcal O}({h^2},\Delta {t^2}),$$

(15)

where overdot denotes differentiation with respect to time.

For the configuration of particles shown in Figure 1, the relative frequency shifts y₁₂ and y₂₁ given by Eqs. (12) can be written, by virtue of the formula (15), in the form

$${y_{12}}(t) = {y_{21}}(t) = - {{{L_3}} \over {2c}}{\rm{n}}_3^{\mathrm{T}} \cdot {\rm{\dot H}}\left({t + {{{{\rm{k}}^{\rm{T}}} \cdot {{\rm{x}}_1}} \over c}} \right) \cdot {{\rm{n}}_3} + \mathcal{O}({h^2},\Delta {t^2}),$$

(16)

so that they are equal to each other up to terms ${\mathcal O({h^2},\Delta {t^2})}$. The photon’s total round-trip frequency shift y₁₂₁ [cf. Eq. (13)] is thus equal to

$${y_{121}}(t) = - {{{L_3}} \over {2c}}{\rm{n}}_3^{\rm{T}}\cdot{\rm{\dot H}}\left({t + {{{{\rm{k}}^{\rm{T}}}\cdot{{\rm{x}}_1}} \over c}} \right)\cdot{{\rm{n}}_3} + {\mathcal O}({h^2},\Delta {t^2}).$$

(17)

There are important cases where the long-wavelength approximation is not valid. These include satellite Doppler tracking measurements and the space-borne LISA detector for gravitational-wave frequencies larger than a few mHz.

2.3 Solar-system-based detectors

Real gravitational-wave detectors do not stay at rest with respect to the TT coordinate system related to the passing gravitational wave, because they also move in the gravitational field of the solar system bodies, as in the case of the LISA spacecraft, or are fixed to the surface of the Earth, as in the case of Earth-based laser interferometers or resonant bar detectors. Let us choose the origin O of the TT coordinate system employed in Section 2.1 to coincide with the solar system barycenter (SSB). The motion of the detector with respect to the SSB will modulate the gravitational-wave signal registered by the detector. One can show that as far as the velocities of the particles (modeling the detector’s parts) with respect to the SSB are non-relativistic, which is the case for all existing or planned detectors, Eqs. (12) can still be used, provided the 3-vectors x_a and n_a (a = 1, 2, 3) will be interpreted as made of the time-dependent components describing the motion of the particles with respect to the SSB.

It is often convenient to introduce the proper reference frame of the detector with coordinates $({{\hat x}^a})$. Because the motion of this frame with respect to the SSB is non-relativistic, we can assume that the transformation between the SSB-related coordinates (x^α) and the detector’s proper reference frame coordinates $({{\hat x}^a})$ has the form

$$\hat t = t,\quad\quad{\hat x^i}(t,{x^k}) = \hat x_{\hat O}^i(t) + O_j^i(t){x^j},$$

(18)

where the functions ${{\hat x}^i}_{\hat O}(t)$ describe the motion of the origin Ô of the proper reference frame with respect to the SSB, and the functions ${\rm{O}}_j^i(t)$ account for the different orientations of the spatial axes of the two reference frames. One can compute some of the quantities entering Eqs. (12) in the detector’s coordinates rather than in the TT coordinates. For instance, the matrix Ĥ of the wave-induced spatial metric perturbation in the detector’s coordinates is related to the matrix H of the spatial metric perturbation produced by the wave in the TT coordinate system through equation

$${\rm{\hat H}}(t) = {({\rm{O}}{(t)^{- 1}})^{\rm{T}}} \cdot {\rm{H}}(t) \cdot {\rm{O}}{(t)^{- 1}},$$

(19)

where the matrix O has elements ${\rm{O}}_j^i$. If the transformation matrix O is orthogonal, then O⁻¹ = OT, and Eq. (19) simplifies to

$${\rm{\hat H}}(t) = {\rm{O}}(t) \cdot {\rm{H}}(t) \cdot {\rm{O}}{(t)^{\rm{T}}}.$$

(20)

See [31, 54, 68, 78] for more details.

For a ground-based laser-interferometric detector, the long-wavelength approximation can be employed (however, see [27, 115, 114] for a discussion of importance of high-frequency corrections, which modify the interferometer response function computed within the long-wavelength approximation). In the case of an interferometer in a standard Michelson and equal-arm configuration (such configurations can be represented by Figure 1 with the particle 1 corresponding to the corner station of the interferometer and with L₂ = L₃ = L), the observed relative frequency shift Δν(t)/ν₀ is equal to the difference of the round-trip frequency shifts in the two detector’s arms [136]:

$${{\Delta \nu (t)} \over {{\nu _0}}} = {y_{131}}(t) - {y_{121}}(t).$$

(21)

Let (x_d, y_d, z_d) be the components (with respect to the TT coordinate system) of the 3-vector r_d connecting the origin of the TT coordinate system with the corner station of the interferometer. Then x₁ = (x_d, y_d, z_d)^T, k^T · x₁ = −z_d, and, making use of Eq. (17), the relative frequency shift (21) can be written as

$${{\Delta \nu (t)} \over {{\nu _0}}} = {L \over c}\left({{\rm{n}}_2^T \cdot {\rm{\dot H}}\left({t - {{{z_d}} \over c}} \right) \cdot {{\rm{n}}_2} - {\rm{n}}_3^{\rm{T}} \cdot {\rm{\dot H}}\left({t - {{{z_d}} \over c}} \right) \cdot {{\rm{n}}_3}} \right).$$

(22)

The difference Δϕ(t) of the phase fluctuations measured, say, by a photo detector, is related to the corresponding relative frequency fluctuations Δν(t) by

$${{\Delta \nu (t)} \over {{\nu _0}}} = {1 \over {2\pi {\nu _0}}}{{{\rm{d}}\Delta \phi (t)} \over {{\rm{d}}t}}.$$

(23)

One can integrate Eq. (22) to write the phase change Δϕ(t) as

$$\Delta \phi (t) = 4\pi {\nu _0}Lh(t),$$

(24)

where the dimensionless function h,

$$h(t):={1 \over 2}\left({{\rm{n}}_2^{\rm{T}} \cdot {\rm{H}}\left({t - {{{z_{\rm{d}}}} \over c}} \right) \cdot {{\rm{n}}_2} - {\rm{n}}_3^{\rm{T}} \cdot {\rm{H}}\left({t - {{{z_{\rm{d}}}} \over c}} \right) \cdot {{\rm{n}}_3}} \right),$$

(25)

is the response function of the interferometric detector to a plane gravitational wave in the long-wavelength approximation. To get Eqs. (24)–(25) directly from Eqs. (22)–(23) one should assume that the quantities n₂, n₃, and z_d [entering Eq. (22)] do not depend on time t. But the formulae (24)–(25) can also be used in the case when those quantities are time dependent, provided the velocities ẋ_a of the detector’s parts with respect to the SSB are non-relativistic. The error we make in such cases is on the order of ${\mathcal O(hv)}$, where v is a typical value of the velocities ẋ_a. Thus, the response function of the Earth-based interferometric detector equals

$$h(t) = {1 \over 2}\left({{{\rm{n}}_2}{{(t)}^{\rm{T}}}\cdot{\rm{H}}\left({t - {{{z_{\rm{d}}}(t)} \over c}} \right)\cdot{{\rm{n}}_2}(t) - {{\rm{n}}_3}{{(t)}^{\rm{T}}}\cdot{\rm{H}}\left({t - {{{z_{\rm{d}}}(t)} \over c}} \right)\cdot{{\rm{n}}_3}(t)} \right),$$

(26)

where all quantities here are computed in the SSB-related TT coordinate system. The same response function can be written in terms of a detector’s proper-reference-frame quantities as follows

$$h(t) = {1 \over 2}\left({{\rm{\hat n}}_2^{\rm{T}}\cdot{\rm{\hat H}}\left({t - {{{z_{\rm{d}}}(t)} \over c}} \right)\cdot{{{\rm{\hat n}}}_2} - {\rm{\hat n}}_3^{\rm{T}}\cdot{\rm{\hat H}}\left({t - {{{z_{\rm{d}}}(t)} \over c}} \right)\cdot{{{\rm{\hat n}}}_3}} \right),$$

(27)

where the matrices Ĥ and H are related to each other by means of formula (20). In Eq. (27) the proper-reference-frame components ${{{\rm{\hat n}}}_2}$ and ${{{\rm{\hat n}}}_3}$ of the unit vectors directed along the interferometer arms can be treated as constant (i.e., time independent) quantities.

From Eqs. (26) and (3) it follows that the response function h is a linear combination of the two wave polarizations h₊ and h_×, so it can be written as

$$h(t) = {F_ +}(t){h_ +}\left({t - {{{z_{\rm{d}}}(t)} \over c}} \right) + {F_ \times}(t){h_ \times}\left({t - {{{z_{\rm{d}}}(t)} \over c}} \right).$$

(28)

The functions F₊ and F_× are the interferometric beam-pattern functions. They depend on the location of the detector on Earth, the position of the gravitational-wave source in the sky, and the polarization angle of the wave (this angle describes the orientation, with respect to the detector, of the axes relative to which the plus and cross polarizations of the wave are defined, see, e.g., Figure 9.2 in [135]). Derivation of the explicit formulae for the interferometric beam patterns F₊ and F_× can be found, e.g., in Appendix C of [66].

In the long-wavelength approximation, the response function of the interferometric detector can be derived directly from the equation of geodesic deviation [126]. Then the response is defined as the relative difference between the wave-induced changes of the proper lengths of the two arms, i.e., $h(t): = (\Delta {{\hat L}_2}(t) - \Delta {{\hat L}_3}(t))/{{\hat L}_0}$, where ${{\hat L}_0} + \Delta {{\hat L}_2}(t)$ and ${{\hat L}_0} + \Delta {{\hat L}_3}(t)$ are the instantaneous values of the proper lengths of the interferometer’s arms and ${{\hat L}_0}$ is the unperturbed proper length of these arms.

In the case of an Earth-based resonant-bar detector the long-wavelength approximation is very accurate and the dimensionless detector’s response function can be defined as $h(t): = \Delta \hat L(t)/{{\hat L}_0}$, where $\Delta \hat L(t)$ is the wave-induced change of the proper length ${{\hat L}_0}$ of the bar. The response function h can be computed in terms of the detector’s proper-reference-frame quantities from the formula (see, e.g., Section 9.5.2 in [135])

$$h(t) = {{\rm{\hat n}}^T} \cdot {\rm{\hat H}}\left({t - {{{z_{\rm{d}}}(t)} \over c}} \right) \cdot {\rm{\hat n}},$$

(29)

where the column matrix ${{\rm{\hat n}}}$ consists of the components (computed in the proper reference frame of the detector) of the unit vector n directed along the symmetry axis of the bar. The response function (29) can be written as a linear combination of the wave polarizations h₊ and h_×, i.e., the formula (28) is also valid for the resonant-bar response function but with some bar-pattern functions F₊ and F_× different from the interferometric beam-pattern functions. Derivation of the explicit form of the bar patterns can be found, e.g., in Appendix C of [66].

2.4 General form of the response function

In many cases of interest the response function of the detector to a plane gravitational wave can be written as a linear combination of n waveforms h_k(t; ξ) which all depend on the same set of parameters ξ = (ξ₁, …, ξ_m)] with constant amplitudes a_k(k = 1, …, n),

$$h(t;\theta) = \sum\limits_{k = 1}^n {{a_k}{h_k}(t;\xi)},$$

(30)

where the vector θ collects all the signal’s parameters. It is convenient to introduce column matrices

$${\rm a}: = \left({\begin{array}{*{20}c} {{a_1}} \\ \vdots \\ {{a_n}} \\ \end{array}} \right),\, {\rm h}(t;\xi): = \left({\begin{array}{*{20}c} {{h_1}(t;\xi)} \\ \vdots \\ {{h_n}(t;\xi)} \\ \end{array}} \right).$$

(31)

then one can briefly write θ = (a, ξ) and the response (30) can be written as

$$h(t;\theta) = {{\rm{a}}^{\rm{T}}} \cdot {\rm{h}}(t;\xi).$$

(32)

The functions h_k (k = 1, …, n) are independent of the parameters a. The parameters a are called extrinsic (or amplitude) parameters whereas the parameters ξ are called intrinsic.

Eq. (32) with n = 4 is a model of the response of the space-based detector LISA to gravitational waves from a binary system [78]. Also for n = 4 the same equation models the response of a ground-based detector to a continuous source of gravitational waves like a rotating neutron star [68]. For ground-based detectors the long-wavelength approximation can be applied and within this approximation the functions h_k (k = 1, …, 4) are given by

$$\begin{array}{*{20}c} {{h_1}(t;\xi) = u(t;\xi)\cos \phi (t;\xi),} \\ {{h_2}(t;\xi) = v(t;\xi)\cos \phi (t;\xi),} \\ {{h_3}(t;\xi) = u(t;\xi)\sin \phi (t;\xi),} \\ {{h_4}(t;\xi) = v(t;\xi)\sin \phi (t;\xi),} \\ \end{array}$$

(33)

where ϕ(t; ξ) is the phase modulation of the signal and u(t; ξ) and v(t; ξ) are slowly varying amplitude modulations. The gravitational-wave signal from spinning neutron stars may consist of several components of the form (32). For short observation times over which the amplitude modulation functions are nearly constant, the response of the ground-based detector can further be approximated by

$${R_j}(\delta) = {C_{0j}}{P_j}(R) + {C_{1j}}{P_j}({R{\prime}}),\,j = 0,1,$$

(34)

where A₀ and ϕ₀ are constant amplitude and initial phase, respectively, and g(t; ξ) is a slowly varying function of time. Eq. (34) is a good model for the response of a detector to the gravitational wave from a coalescing compact binary system [135, 30]. We would like to stress that not all deterministic gravitational-wave signals may be cast into the general form (32).

3 Statistical Theory of Signal Detection

The gravitational-wave signal will be buried in the noise of the detector and the data from the detector will be a random (or stochastic) process. Consequently, the problem of extracting the signal from the noise is a statistical one. The basic idea behind signal detection is that the presence of the signal changes the statistical characteristics of the data x, in particular its probability distribution. When the signal is absent the data have probability density function (pdf) p₀(x), and when the signal is present the pdf is p₁(x).

A thorough introduction to probability theory and mathematical statistics can be found, e.g., in [51, 130, 131, 101]. A full exposition of statistical theory of signal detection that is only outlined here can be found in the monographs [147, 74, 143, 139, 87, 58, 107]. A general introduction to stochastic processes is given in [145] and advanced treatment of the subject can be found in [84, 146]. A concise introduction to the statistical theory of signal detection and time series analysis is contained in Chapters 3 and 4 of [66].

3.1 Hypothesis testing

The problem of detecting the signal in noise can be posed as a statistical hypothesis testing problem. The null hypothesis H₀ is that the signal is absent from the data and the alternative hypothesis H₁ is that the signal is present. A hypothesis test (or decision rule) δ is a partition of the observation set into two subsets, ${\mathcal R}$ and its complement ${\mathcal R\prime}$. If data are in ${\mathcal R}$ we accept the null hypothesis, otherwise we reject it. There are two kinds of errors that we can make. A type I error is choosing hypothesis H₁ when H₀ is true and a type II error is choosing H₀ when H₁ is true. In signal detection theory the probability of a type I error is called the false alarm probability, whereas the probability of a type II error is called the false dismissal probability. 1 − (false dismissal probability) is the probability of detection of the signal. In hypothesis testing theory, the probability of a type I error is called the significance of the test, whereas 1 − (probability of type II error) is called the power of the test.

The problem is to find a test that is in some way optimal. There are several approaches to finding such a test. The subject is covered in detail in many books on statistics, for example, see [72, 51, 80, 83].

3.1.1 Bayesian approach

In the Bayesian approach we assign costs to our decisions; in particular we introduce positive numbers C_ij, i, j = 0, 1, where C_ij is the cost incurred by choosing hypothesis H_i when hypothesis H_j is true. We define the conditional risk R of a decision rule δ for each hypothesis as

$${R_j}(\delta) = {C_{0j}}{P_j}(\mathcal{R}) + {C_{1j}}{P_j}({\mathcal{R}{\prime}}), \quad\quad j = 0,1,$$

(35)

where P_j is the probability distribution of the data when hypothesis H_j is true. Next, we assign probabilities π₀ and π₁ = 1 − π₀ to the occurrences of hypotheses H₀ and H₁, respectively. These probabilities are called a priori probabilities or priors. We define the Bayes risk as the overall average cost incurred by the decision rule δ:

$$r(\delta) = {\pi _0}{R_0}(\delta) + {\pi _1}{R_1}(\delta).$$

(36)

Finally we define the Bayes rule as the rule that minimizes the Bayes risk r(δ).

3.1.2 Minimax approach

Very often in practice we do not have control over or access to the mechanism generating the state of nature and we are not able to assign priors to various hypotheses. In such a case one criterion is to seek a decision rule that minimizes, over all δ, the maximum of the conditional risks, R₀(δ) and R₁(δ). A decision rule that fulfills that criterion is called a minimax rule.

3.1.3 Neyman-Pearson approach

In many problems of practical interest the imposition of a specific cost structure on the decisions made is not possible or desirable. The Neyman-Pearson approach involves a trade-off between the two types of errors that one can make in choosing a particular hypothesis. The Neyman-Pearson design criterion is to maximize the power of the test (probability of detection) subject to a chosen significance of the test (false alarm probability).

3.1.4 Likelihood ratio test

It is remarkable that all three very different approaches — Bayesian, minimax, and Neyman-Pearson — lead to the same test called the likelihood ratio test [44]. The likelihood ratio Λ is the ratio of the pdf when the signal is present to the pdf when it is absent:

$$\Lambda (x): = {{{p_1}(x)} \over {{p_0}(x)}}.$$

(37)

We accept the hypothesis H₁ if Λ > k, where k is the threshold that is calculated from the costs C_ij, priors π_i, or the significance of the test depending on which approach is being used.

3.2 The matched filter in Gaussian noise

Let h be the gravitational-wave signal we are looking for and let n be the detector’s noise. For convenience we assume that the signal h is a continuous function of time t and that the noise n is a continuous random process. Results for the discrete-in-time data that we have in practice can then be obtained by a suitable sampling of the continuous-in-time expressions. Assuming that the noise is additive the dat x can be written as

$$x(t) = n(t) + h(t).$$

(38)

The autocorrelation function of the noise n is defined as

$${K_n}(t,t{\prime}): = {\rm E}[n(t)n({t{\prime}})],$$

(39)

where E denotes the expectation value.

Let us further assume that the detector’s noise n is a zero-mean and Gaussian random process. It can then be shown that the logarithm of the likelihood function is given by the following Cameron-Martin formula

$$\log \Lambda [x] = \int\nolimits_0^{{T_o}} {q(t)x(t){\rm{d}}t - {1 \over 2}\int\nolimits_0^{{T_o}} {q(t)h(t){\rm{d}}t,}}$$

(40)

where 〈0; T_o〉 is the time interval over which the data was collected and the function q is the solution of the integral equation

$$h(t) = \int\nolimits_0^{{T_o}} {{K_n}(t,{t\prime})q({t\prime}){\rm{d}}t.}$$

(41)

For stationary noise, its autocorrelation function (39) depends on times t and t^′ only through the difference t − t^′. It implies that there exists then an even function κ_n of one variable such that

$${\rm E}[n(t)n({t{\prime}})] = {\kappa _n}(t - {t{\prime}}).$$

(42)

Spectral properties of stationary noise are described by its one-sided spectral density, defined for non-negative frequencies f ≥ 0 through relation

$${S_n}(f) = 2\int\nolimits_{- \infty}^\infty {{\kappa _n}(t){e^{2\pi i ft}}{\rm{d}}t.}$$

(43)

for negative frequencies f < 0, by definition, S_n(f) = 0. The spectral density S_n can also be determined by correlations between the Fourier components of the detector’s noise

$${\rm{E}}[\tilde n(f){\tilde n^{\ast}}({f{\prime}})] = {1 \over 2}{S_n}(\vert f\vert)\delta (f - {f{\prime}}),\, - \infty < f,{f{\prime}} < \infty.$$

(44)

For the case of stationary noise with one-sided spectral density S_n, it is convenient to define the scalar product (x∣y) between any two waveforms x and y,

$$(x\vert y): = 4\Re \int\nolimits_0^\infty {{{\tilde x(f){{\tilde y}^ {\ast}}(f)} \over {{S_n}(f)}}{\rm{d}}f,}$$

(45)

where ℜ denotes the real part of a complex expression, the tilde denotes the Fourier transform, and the asterisk is complex conjugation. Making use of this scalar product, the log likelihood function (40) can be written as

$$\log \Lambda [x] = (x\vert h) - {1 \over 2}(h\vert h).$$

(46)

From the expression (46) we see immediately that the likelihood ratio test consists of correlating the data x with the signal h that is present in the noise and comparing the correlation to a threshold. Such a correlation is called the matched filter. The matched filter is a linear operation on the data.

An important quantity is the optimal signal-to-noise ratio ρ defined by

$$\rho: = \sqrt {(h\vert h).}$$

(47)

By means of Eq. (45) it can be written as

$${\rho ^2} = 4\int\nolimits_0^\infty {{{\vert \tilde h(f){\vert ^2}} \over {{S_n}(f)}}{\rm{d}}f.}$$

(48)

We see in the following that ρ determines the probability of detection of the signal. The higher the signal-to-noise ratio the higher the probability of detection.

An interesting property of the matched filter is that it maximizes the signal-to-noise ratio over all linear filters [44]. This property is independent of the probability distribution of the noise.

3.3 Parameter estimation

Very often we know the waveform of the signal that we are searching for in the data in terms of a finite number of unknown parameters. We would like to find optimal procedures of estimating these parameters. An estimator of a parameter θ is a function $\hat \theta (x)$ that assigns to each data set the “best” guess of the true value of θ. Note that because $\hat \theta (x)$ depends on the random data it is a random variable. Ideally we would like our estimator to be (i) unbiased, i.e., its expectation value to be equal to the true value of the parameter, and (ii) of minimum variance. Such estimators are rare and in general difficult to find. As in the signal detection there are several approaches to the parameter estimation problem. The subject is exposed in detail in [81, 82]. See also [149] for a concise account.

3.3.1 Bayesian estimation

We assign a cost function C(θ^′, θ) of estimating the true value of θ as θ^′. We then associate with an estimator ${\hat \theta}$ a conditional risk or cost averaged over all realizations of data x for each value of the parameter θ:

$${R_\theta}(\hat \theta) = {{\rm{E}}_\theta}[C(\hat \theta, \theta)] = \int\nolimits_X {C(\hat \theta (x),\theta)p(x,\theta)} {\rm{d}}x,$$

(49)

where X is the set of observations and p(x, θ) is the joint probability distribution of data x and parameter θ. We further assume that there is a certain a priori probability distribution π(θ) of the parameter θ. We then define the Bayes estimator as the estimator that minimizes the average risk defined as

$$r(\hat \theta) = {\rm{E}}[{R_\theta}(\hat \theta)] = \int\nolimits_X {\int\nolimits_\Theta {C(\hat \theta (x),\theta)} p(x,\theta)} \pi (\theta){\rm{d}}\theta {\rm{d}}x,$$

(50)

where E is the expectation value with respect to an a priori distribution π, and Θ is the set of observations of the parameter θ. It is not difficult to show that for a commonly used cost function

$$C({\theta {\prime}},\theta) = {({\theta {\prime}},\theta)^2},$$

(51)

the Bayesian estimator is the conditional mean of the parameter θ given data x, i.e.,

$$\hat \theta (x) = E[\theta \vert x] = \int\nolimits_\Theta {\theta p(\theta \vert x)d\theta,}$$

(52)

where p(θ∣x) is the conditional probability density of parameter θ given the data x.

3.3.2 Maximum a posteriori probability estimation

Suppose that in a given estimation problem we are not able to assign a particular cost function C(θ′, θ). Then a natural choice is a uniform cost function equal to 0 over a certain interval I_θ of the parameter θ. From Bayes theorem [28] we have

$$p(\theta \vert x) = {{p(x,\theta)\pi (\theta)} \over {p(x)}},$$

(53)

where p(x) is the probability distribution of data x. Then, from Eq. (50) one can deduce that for each data point x the Bayes estimate is any value of θ that maximizes the conditional probability p(θ∣x). The density p(θ∣x) is also called the a posteriori probability density of parameter θ and the estimator that maximizes p(θ∣x) is called the maximum a posteriori (MAP) estimator. It is denoted by ${{\hat \theta}_{{\rm{MAP}}}}$. We find that the MAP estimators are solutions of the following equation

$${{\partial \log p(x,\theta)} \over {\partial \theta}} = - {{\partial \log \pi (\theta)} \over {\partial \theta}},$$

(54)

which is called the MAP equation.

3.3.3 Maximum likelihood estimation

Often we do not know the a priori probability density of a given parameter and we simply assign to it a uniform probability. In such a case maximization of the a posteriori probability is equivalent to maximization of the probability density p(x, θ) treated as a function of θ. We call the function l(θ, x) ≔ p(x, θ) the likelihood function and the value of the parameter θ that maximizes l(θ, x) the maximum likelihood (ML) estimator. Instead of the function l we can use the function Λ(θ, x) = l(θ, x)/p(x) (assuming that p(x) > 0). Λ is then equivalent to the likelihood ratio [see Eq. (37)] when the parameters of the signal are known. Then the ML estimators are obtained by solving the equation

$${{\partial \log \Lambda (\theta, x)} \over {\partial \theta}} = 0,$$

(55)

which is called the ML equation.

3.4 Fisher information and Cramèr-Rao bound

It is important to know how good our estimators are. We would like our estimator to have as small a variance as possible. There is a useful lower bound on variances of the parameter estimators called the Cramèr-Rao bound, which is expressed in terms of the Fisher information matrix Γ(θ) For the signal h(t; θ), which depends on K parameters collected into the vector θ = (θ₁,…, θ_k), the components of the matrix Γ(θ) are defined as

$$\Gamma {(\theta)_{ij}}: = {\rm{E}}\left[ {{{\partial \log \Lambda [x;\theta ]} \over {\partial {\theta _i}}}{{\partial \log \Lambda [x;\theta ]} \over {\partial {\theta _j}}}} \right] = - {\rm{E}}\left[ {{{{\partial ^2}\log \Lambda [x;\theta ]} \over {\partial {\theta _i}\partial {\theta _j}}}} \right],\,i,j = 1, \ldots, K.$$

(56)

The Cramèr-Rao bound states that for unbiased estimators the covariance matrix C(θ) of the estimators θ fulfills the inequality

$${\rm{C}}(\theta) \geq\Gamma {(\theta)^{- 1}}.$$

(57)

(The inequality A ≥ B for matrices means that the matrix A − B is nonnegative definite.)

A very important property of the ML estimators is that asymptotically (i.e., for a signal-to-noise ratio tending to infinity) they are (i) unbiased, and (ii) they have a Gaussian distribution with covariance matrix equal to the inverse of the Fisher information matrix.

In the case of Gaussian noise the components of the Fisher matrix are given by

$$\Gamma {(\theta)_{ij}} = \left({\left. {{{\partial h(t;\theta)} \over {\partial {\theta _i}}}} \right\vert {{\partial h(t;\theta)} \over {\partial {\theta _j}}}} \right),\,i,j = 1, \ldots, K,$$

(58)

where the scalar product (·∣·) is defined in Eq. (45).

4 Maximum-likelihood Detection in Gaussian Noise

In this section, we study the detection of a deterministic gravitational-wave signal h(t; θ) of the general form given by Eq. (32) and the estimation of its parameters θ using the maximum-likelihood (ML) principle. We assume that the noise n(t) in the detector is a zero-mean, Gaussian, and stationary random process. The data x in the detector, in the case when the gravitational-wave signal h(t; θ) is present, is x(t; θ) = n(t) + h(t; θ). The parameters θ = (a, ξ) of the signal (32) split into extrinsic (or amplitude) parameters a and intrinsic ones ξ.

4.1 The ${\mathcal F}$-statistic

For the gravitational-wave signal h(t; a, ξ) of the form given in Eq. (32) the log likelihood function (46) can be written as

$$\log \Lambda [x;{\rm{a}},\xi ] = {{\rm{a}}^{\rm{T}}} \cdot {\rm{N}}[x;\xi ] - {1 \over 2}{{\rm{a}}^{\rm{T}}} \cdot {\rm{M(}}\xi) \cdot {\rm{a,}}$$

(59)

where the components of the column n × 1 matrix N and the square n × n matrix M are given by

$${N_k}[x;\xi ]: = (x\vert {h_k}(t;\xi)),\,{M_{kl}}(\xi): = ({h_k}(t;\xi)\vert {h_l}(t;\xi)),\,k,l = 1, \ldots, n.$$

(60)

The ML equations for the extrinsic parameters a, ∂ logΛ[x; a, ξ]/∂a = 0, can be solved explicitly to show that the ML estimators â of the parameters a are given by

$${\rm{\hat a}}[x;\xi ] = {\rm{M}}{(\xi)^{- 1}} \cdot {\rm{N}}[x;\xi ].$$

(61)

Replacing the extrinsic parameters a in Eq. (59) by their ML estimators â, we obtain the reduced log likelihood function,

$${\mathcal F}[x;\xi ]: = \log \Lambda [x;{\rm{\hat a}}[x;\xi ],\xi ] = {1 \over 2}{\rm N}{[x;\xi ]^{\rm{T}}} \cdot {\rm{M}}{(\xi)^{- 1}} \cdot {\rm{N}}[x;\xi ],$$

(62)

that we call the ${\mathcal F}$-statistic. The ${\mathcal F}$-statistic depends nonlinearly on the intrinsic parameters ξ of the signal, it does not depend on the extrinsic parameters a.

The procedure to detect the gravitational-wave signal of the form (32) and estimate its parameters consists of two parts. The first part is to find the (local) maxima of the ${\mathcal F}$-statistic (62) in the intrinsic parameters space. The ML estimators ${\hat \xi}$ of the intrinsic parameters ξ are those values of ξ for which the ${\mathcal F}$-statistic attains a maximum. The second part is to calculate the estimators â of the extrinsic parameters a from the analytic formula (61), where the matrix M and the correlations N are calculated for the parameters ξ replaced by their ML estimators ${\hat \xi}$ obtained from the first part of the analysis. We call this procedure the maximum likelihood detection. See Section 4.6 for a discussion of the algorithms to find the (local) maxima of the ${\mathcal F}$-statistic.

4.1.1 Targeted searches

The ${\mathcal F}$-statistic can also be used in the case when the intrinsic parameters are known. An example of such an analysis called a targeted search is the search for a gravitational-wave signal from a known pulsar. In this case assuming that gravitational-wave emission follows the radio timing, the phase of the signal is known from pulsar observations and the only unknown parameters of the signal are the amplitude (or extrinsic) parameters a [see Eq. (30)]. To detect the signal one calculates the ${\mathcal F}$-statistic for the known values of the intrinsic parameters and compares it to a threshold [67]. When a statistically-significant signal is detected, one then estimates the amplitude parameters from the analytic formulae (61).

In [109] it was shown that the maximum-likelihood ${\mathcal F}$-statistic can be interpreted as a Bayes factor with a simple, but unphysical, amplitude prior (and an additional unphysical sky-position weighting). Using a more physical prior based on an isotropic probability distribution for the unknown spin-axis orientation of emitting systems, a new detection statistic (called the ${\mathcal B}$-statistic) was obtained. Monte Carlo simulations for signals with random (isotropic) spin-axis orientations show that the ${\mathcal B}$-statistic is more powerful (in terms of its expected detection probability) than the ${\mathcal F}$-statistic. A modified version of the ${\mathcal F}$-statistic that can be more powerful than the original one has been studied in [20].

4.2 Signal-to-noise ratio and the Fisher matrix

The detectability of the signal h(t; θ) is determined by the signal-to-noise ratio ρ. In general it depends on all the signal’s parameters θ and can be computed from [see Eq. (47)]

$$\rho (\theta) = \sqrt {(h(t;\theta)\vert h(t;\theta)}).$$

(63)

The signal-to-noise ratio for the signal (32) can be written as

$$\rho ({\rm{a}},\xi) = \sqrt {{{\rm{a}}^{\rm{T}}} \cdot {\rm{M}}(\xi) \cdot {\rm{a}}},$$

(64)

where the components of the matrix M(ξ) are defined in Eq. (60).

The accuracy of estimation of the signal’s parameters is determined by Fisher information matrix Γ. The components of Γ in the case of the Gaussian noise can be computed from Eq. (58). For the signal given in Eq. (32) the signal’s parameters (collected into the vector θ) split into extrinsic and intrinsic parameters: θ = (a, ξ), where a = (a₁, …, a_n) and ξ = (ξ₁,…, ξ_m). It is convenient to distinguish between extrinsic and intrinsic parameter indices. Therefore, we use calligraphic lettering to denote the intrinsic parameter indices: ${\xi _{\mathcal A}},{\mathcal A} = 1, \ldots,m$. The matrix Γ has dimension (n + m) × (n + m) and it can be written in terms of four block matrices for the two sets of the parameters a and ξ,

$$\Gamma ({\rm a},\xi) = \left({\begin{array}{*{20}c} {{\Gamma _{{\rm{aa}}}}(\xi)} & {{\Gamma _{{\rm a}\xi}}({\rm a},\xi)} \\ {{\Gamma _{{\rm a}\xi}}{{({\rm a},\xi)}^{\rm{T}}}} & {{\Gamma _{\xi \xi}}({\rm a},\xi)} \\ \end{array}} \right),$$

(65)

where Γ_aa is an n × n matrix with components (∂h/∂a_i∣∂h/∂a_j) (i, j = 1, …, n), Γ_aξ is an n × m matrix with components $(\partial h/\partial {a_i}\backslash \partial h/\partial {\xi _\mathcal A})(i = 1, \ldots, n,\mathcal A = 1, \ldots, m)$, and finally Γ_ξξ is m × m matrix with components $(\partial h/\partial {\xi _\mathcal A}\backslash \partial h/\partial {\xi _\mathcal B})(\mathcal A,\mathcal B = 1, \ldots, m)$.

We introduce two families of the auxiliary n × n square matrices ${{\rm{F}}_{(\mathcal A)}}$ and ${{\rm{S}}_{({\mathcal A}{\mathcal B})}}({\mathcal A},{\mathcal B} = 1, \ldots, m)$, which depend on the intrinsic parameters ξ only (the indexes ${\mathcal A},{\mathcal B}$ within parentheses mean that they serve here as the matrix labels). The components of the matrices ${{\rm{F}}_{(\mathcal A)}}$ and ${{\rm{S}}_{({\mathcal A}{\mathcal B})}}$ are defined as follows:

$${F_{{{(\mathcal{A})}^{ij}}}}(\xi): = \left({{h_i}(t;\xi)\left\vert {{{\partial {h_j}(t;\xi)} \over {\partial {\xi _\mathcal{A}}}}} \right.} \right),\,i,j = 1, \ldots, n, \,\mathcal{A} = 1, \ldots, m.$$

(66)

$${S_{{{({\mathcal A}{\mathcal B})}^{ij}}}}(\xi): = \left({{{\partial {h_i}(t;\xi)} \over {\partial {\xi _{\mathcal A}}}}\left\vert {{{\partial {h_j}(t;\xi)} \over {\partial {\xi _{\mathcal B}}}}} \right.} \right),\,i,j = 1, \ldots, n,\,{\mathcal A},{\mathcal B} = 1, \ldots, m.$$

(67)

Making use of the definitions (60) and (66)–(67) one can write the more explicit form of the matrices Γ_aa, Γ_aξ, and Γ_ξξ,

$${\Gamma _{aa}}(\xi) = {\rm{M}}(\xi),$$

(68)

$${\Gamma _{a\xi}}({\rm{a}},\xi) = ({{\rm{F}}_{(1)}}(\xi) \cdot {\rm{a}} \cdots {{\rm{F}}_{(m)}}(\xi) \cdot {\rm{a}}),$$

(69)

$${\Gamma _{\xi \xi}}({\rm{a}},\xi) = \left({\begin{array}{*{20}c} {{{\rm{a}}^{\rm{T}}} \cdot {{\rm{S}}_{(11)}}(\xi) \cdot {\rm{a}} \cdots {{\rm{a}}^{\rm{T}}} \cdot {{\rm{S}}_{(1m)}}(\xi) \cdot {\rm{a}}} \\ {\cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots} \\ {{{\rm{a}}^{\rm{T}}} \cdot {{\rm{S}}_{(m1)}}(\xi) \cdot {\rm{a}} \cdots {{\rm{a}}^{\rm{T}}} \cdot {{\rm{S}}_{(mm)}}(\xi) \cdot {\rm{a}}} \\ \end{array}} \right).$$

(70)

The notation introduced above means that the matrix Γ_aξ can be thought of as a 1 × m row matrix made of n × 1 column matrices ${{\rm{F}}_{(\mathcal A)}}$. Thus, the general formula for the component of the matrix Γ_aξ is

$${({\Gamma _{{\rm{a}}\xi}})_{i\mathcal{A}}} = ({{\rm{F}}_{(\mathcal{A)}} \cdot {\rm{a}})_i} = \sum\limits_{j = 1}^n {{F_{(\mathcal{A})ij}}{a_j}, \quad\quad \mathcal{A} = 1, \ldots, m, \quad i = 1, \ldots, n.}$$

(71)

The general component of the matrix Γ_ξξ is given by

$${({\Gamma _{\xi \xi}})_{{\mathcal A}{\mathcal B}}} = {{\rm a}^{\rm T}} \cdot {S_{(\mathcal{A}\mathcal{B})}} \cdot {\rm a} = \sum\limits_{i = 1}^n {\sum\limits_{j = 1}^n {{S_{({\mathcal A}{\mathcal B})ij}}{a_i}{a_j},\,{\mathcal A},\mathcal{B} = 1, \ldots, m.}}$$

(72)

The covariance matrix C, which approximates the expected covariances of the ML estimators of the parameters θ, is defined as Γ⁻¹ Applying the standard formula for the inverse of a block matrix [90] to Eq. (65), one gets

$${\rm C}({\rm{a}},\xi) = \left({\begin{array}{*{20}c} {{{\rm C}_{{\rm{aa}}}}({\rm{a}},\xi)} & {{{\rm C}_{{\rm{a}}\xi}}({\rm{a}},\xi)} \\ {{{\rm C}_{{\rm{a}}\xi}}{{({\rm{a}},\xi)}^T}} & {{{\rm C}_{\xi \xi}}({\rm{a}},\xi)} \\ \end{array}} \right),$$

(73)

where the matrices C_aa, C_aξ, and C_ξξ can be expressed in terms of the matrices Γ_aa = M, Γ_aξ, and Γ_ξξ as follows:

$${{\rm{C}}_{{\rm{aa}}}}({\rm{a}},\xi) = {\rm{M(}}\xi)^{- 1} + {\rm{M(}}\xi)^{- 1} \cdot {\Gamma _{{\rm{a}}\xi}}({\rm{a,}}\xi) \cdot \bar \Gamma {({\rm{a,}}\xi)^{- 1}} \cdot {\Gamma _{{\rm{a}}\xi}}{({\rm{a,}}\xi)^{\rm{T}}} \cdot {\rm{M(}}\xi)^{- 1},$$

(74)

$${{\rm{C}}_{{\rm{a}}\xi}}({\rm{a}},\xi) = - {\rm{M}}{(\xi)^{- 1}}\cdot{\Gamma _{{\rm{a}}\xi}}({\rm{a}},\xi)\cdot\bar \Gamma {({\rm{a}},\xi)^{- 1}},$$

(75)

$${{\rm{C}}_{\xi \xi}}({\rm{a}},\xi) = \bar \Gamma {({\rm{a}},\xi)^{- 1}}.$$

(76)

In Eqs. (74)–(76) we have introduced the m × m matrix:

$$\bar \Gamma ({\rm{a,}}\xi): = {\Gamma _{\xi \xi}}({\rm{a,}}\xi) - {\Gamma _{{\rm{a}}\xi}}{({\rm{a}},\xi)^{\rm{T}}} \cdot {\rm{M(}}\xi)^{- 1} \cdot {\Gamma _{{\rm{a}}\xi}}({\rm{a,}}\xi).$$

(77)

We call the matrix ${\bar \Gamma}$ (which is the Schur complement of the matrix M) the projected Fisher matrix (onto the space of intrinsic parameters). Because the matrix ${\bar \Gamma}$ is the inverse of the intrinsic-parameter submatrix C_ξξ of the covariance matrix C, it expresses the information available about the intrinsic parameters that takes into account the correlations with the extrinsic parameters. The matrix ${\bar \Gamma}$ is still a function of the putative extrinsic parameters.

We next define the normalized projected Fisher matrix (which is the m × m square matrix)

$${\bar \Gamma _n}({\rm{a}},\xi): = {{\bar \Gamma ({\rm{a,}}\xi)} \over {\rho {{({\rm{a,}}\xi)}^2}}},$$

(78)

where ρ is the signal-to-noise ratio. Making use of the definition (77) and Eqs. (71)–(72) we can show that the components of this matrix can be written in the form

$$({\bar \Gamma _n}{({\rm{a}},\xi)_{\mathcal{A}\mathcal{B}}}) = {{{{\rm{a}}^{\rm{T}}} \cdot {{\rm{A}}_{(\mathcal{A}\mathcal{B})}}(\xi) \cdot {\rm{a}}} \over {{{\rm{a}}^{\rm{T}}} \cdot {\rm{M(}}\xi) \cdot {\rm{a}}}},\mathcal{A},\mathcal{B} = 1, \ldots, m,$$

(79)

where ${{\rm{A}}_{({\mathcal A}{\mathcal B})}}$ is the n × n matrix defined as

$${{\rm{A}}_{(\mathcal{A}\mathcal{B})}}(\xi): = {S_{(\mathcal{A}\mathcal{B})}}(\xi) - {{\rm{F}}_{(\mathcal{A})}}{(\xi)^{\rm{T}}} \cdot {\rm{M(}}\xi)^{- 1} \cdot {{\rm{F}}_{(B)}}(\xi),\,\mathcal{A},\mathcal{B} = 1, \ldots, m.$$

(80)

From the Rayleigh principle [90] it follows that the minimum value of the component ${({{\bar \Gamma}_n}({\rm{a}},\xi))_{{\mathcal A}{\mathcal B}}}$ is given by the smallest eigenvalue of the matrix ${{\rm{M}}^{- 1}} \cdot {{\rm{A}}_{({\mathcal A}{\mathcal B})}}$. Similarly, the maximum value of the component ${({{\bar \Gamma}_n}({\rm{a}},\xi))_{{\mathcal A}{\mathcal B}}}$ is given by the largest eigenvalue of that matrix.

Because the trace of a matrix is equal to the sum of its eigenvalues, the m × m square matrix ${\tilde \Gamma}$ with components

$${(\tilde \Gamma (\xi))_{\mathcal{A}\mathcal{B}}}: = {1 \over n}{\rm{Tr(M(}}\xi)^{- 1} \cdot {{\rm{A}}_{(\mathcal{A}\mathcal{B})}}(\xi)),\,\mathcal{A},\mathcal{B} = 1, \ldots, m.$$

(81)

expresses the information available about the intrinsic parameters, averaged over the possible values of the extrinsic parameters. Note that the factor 1/n is specific to the case of n extrinsic parameters. We shall call ${\tilde \Gamma}$ the reduced Fisher matrix. This matrix is a function of the intrinsic parameters alone. We shall see that the reduced Fisher matrix plays a key role in the signal processing theory that we present here. It is used in the calculation of the threshold for statistically significant detection and in the formula for the number of templates needed to do a given search.

For the case of the signal

$$h(t;{A_0},{\phi _0},\xi) = {A_0}g(t;\xi)\cos (\phi (t;\xi) - {\phi _0}),$$

(82)

the normalized projected Fisher matrix ${{\bar \Gamma}_n}$ is independent of the extrinsic parameters A₀ and ϕ₀, and it is equal to the reduced matrix ${\tilde \Gamma}$ [102]. The components of ${\tilde \Gamma}$ are given by

$${\tilde \Gamma _{\mathcal{A}\mathcal{B}}} = {({\Gamma _0})_{\mathcal{A}\mathcal{B}}} - {{{{({\Gamma _0})}_{{\phi _0}\mathcal{A}}}{{({\Gamma _0})}_{{\phi _0}\mathcal{B}}}} \over {{{({\Gamma _0})}_{{\phi _0}{\phi _0}}}}},$$

(83)

where Γ₀ is the Fisher matrix for the signal g(t; ξ) cos (ϕ(t; ξ) − ϕ₀).

4.3 False alarm and detection probabilities

4.3.1 False alarm and detection probabilities for known intrinsic parameters

We first present the false alarm and detection probabilities when the intrinsic parameters ξ of the signal are known. In this case the ${\mathcal F}$-statistic is a quadratic form of the random variables that are correlations of the data. As we assume that the noise in the data is Gaussian and the correlations are linear functions of the data, ${\mathcal F}$ is a quadratic form of the Gaussian random variables. Consequently the ${\mathcal F}$-statistic has a distribution related to the χ² distribution. One can show (see Section III B in [65]) that for the signal given by Eq. (30), $2{\mathcal F}$ has a χ² distribution with n degrees of freedom when the signal is absent and noncentral χ² distribution with n degrees of freedom and non-centrality parameter equal to the square of the signal-to-noise ratio when the signal is present.

As a result the pdfs p₀ and p₁ of the ${\mathcal F}$-statistic, when the intrinsic parameters are known and when respectively the signal is absent or present in the data, are given by

$${p_0}(\mathcal{F}) = {{{\mathcal{F}^{n/2 - 1}}} \over {(n/2 - 1)!}}\exp (- \mathcal{F}),$$

(84)

$${p_1}(\rho,\mathcal{F}) = {{{{(2\mathcal{F})}^{(n/2 - 1)/2}}} \over {{\rho ^{n/2 - 1}}}}{I_{n/2 - 1}}\left({\rho \sqrt {2\mathcal{F}}} \right)\exp \left({- \mathcal{F} - {1 \over 2}{\rho ^2}} \right),$$

(85)

where n is the number of degrees of freedom of χ² distribution and I_n/2−1 is the modified Bessel function of the first kind and order n/2 − 1. The false alarm probability P_F is the probability that ${\mathcal F}$ exceeds a certain threshold ${{\mathcal F}_0}$ when there is no signal. In our case we have

$${P_{\rm{F}}}({\mathcal{F}_0}): = \int\nolimits_{{\mathcal{F}_0}}^\infty {{p_0}(\mathcal{F})d\mathcal{F} = \exp (- {\mathcal{F}_0})\sum\limits_{k = 0}^{n/2 - 1} {{{{\mathcal{F}_0}^k} \over {k!}}.}}$$

(86)

The probability of detection P_D is the probability that ${\mathcal F}$ exceeds the threshold ${{\mathcal F}_0}$ when a signal is present and the signal-to-noise ratio is equal to ρ:

$${P_{\rm{D}}}(\rho, {\mathcal{F}_0}): = \int\nolimits_{{\mathcal{F}_0}}^\infty {p_1}(\rho, \mathcal{F}){\rm{d\mathcal{F}}}.$$

(87)

The integral in the above formula can be expressed in terms of the generalized Marcum Q-function [132, 58], ${P_{\rm{D}}}(\rho, {{\mathcal F}_0}) = Q(\rho, \sqrt {2{{\mathcal F}_0}})$. We see that when the noise in the detector is Gaussian and the intrinsic parameters are known, the probability of detection of the signal depends on a single quantity: the optimal signal-to-noise ratio ρ.

4.3.2 False alarm probability for unknown intrinsic parameters

Next we return to the case in which the intrinsic parameters ξ are not known. Then the statistic ${\mathcal F}[x;\xi ]$ given by Eq. (62) is a certain multiparameter random process called the random field (see monographs [5, 6] for a comprehensive discussion of random fields). If the vector ξ has one component the random field is simply a random process. For random fields we define the autocovariance function ${\mathcal C}$ just in the same way as we define such a function for a random process:

$$\mathcal{C}(\xi, {\xi {\prime}}): = {{\rm{E}}_0}[\mathcal{F}[x;\xi ]\mathcal{F}[x;{\xi {\prime}}]] - {{\rm{E}}_0}[\mathcal{F}[x;\xi ]]{{\rm{E}}_0}[\mathcal{F}[x;{\xi {\prime}}]],$$

(88)

where ξ and ξ^′ are two values of the intrinsic parameter set, and E₀ is the expectation value when the signal is absent. One can show that for the signal (30) the autocovariance function ${\mathcal C}$ is given by

$${\mathcal C}(\xi, {\xi \prime}) = {1 \over 2}{\rm{Tr}}\left({{\rm{Q(}}\xi, {\xi {\prime}}) \cdot {\rm{M(}}{\xi {\prime}})^{- 1} \cdot {\rm{Q(}}\xi,{\xi {\prime}})^{\rm{T}} \cdot {\rm{M(}}\xi)^{- 1}} \right),$$

(89)

where Q is an n × n matrix with components

$${\rm{Q(}}\xi, {\xi {\prime}})_{ij}: = ({h_i}(t;\xi)\vert {h_j}(t;{\xi {\prime}})),\,i,j = 1, \ldots n.$$

(90)

Obviously Q(ξ, ξ) = M(ξ), therefore ${\mathcal C}(\xi, \xi) = n/2$.

One can estimate the false alarm probability in the following way [68]. The autocovariance function ${\mathcal C}$ tends to zero as the displacement Δξ = ξ^′ − ξ increases (it is maximal for Δξ = 0). Thus we can divide the parameter space into elementary cells such that in each cell the autocovariance function ${\mathcal C}$ is appreciably different from zero. The realizations of the random field within a cell will be correlated (dependent), whereas realizations of the random field within each cell and outside of the cell are almost uncorrelated (independent). Thus, the number of cells covering the parameter space gives an estimate of the number of independent realizations of the random field.

We choose the elementary cell with its origin at the point ξ to be a compact region with boundary defined by the requirement that the autocovariance ${\mathcal C}(\xi, {\xi \prime})$ between the origin ξ and any point ξ^′ at the cell’s boundary equals half of its maximum value, i.e., ${\mathcal C}(\xi, \xi)/2$. Thus, the elementary cell is defined by the inequality

$$\mathcal{C}(\xi, {\xi {\prime}}) \leq {1 \over 2}\mathcal{C}(\xi, \xi) = {n \over 4},$$

(91)

with ξ at the cell’s center and ξ^′ on the cell’s boundary.

To estimate the number of cells we perform the Taylor expansion of the autocovariance function up to the second-order terms:

$$\mathcal{C}(\xi ,\xi \prime ) \cong \frac{n}{2} + \sum\limits_{\mathcal{A} = 1}^m {{{\left. {\frac{{\partial \mathcal{C}(\xi ,\xi \prime )}}{{\partial {\xi _\mathcal{A}}\prime }}} \right|}_{\xi \prime= \xi }}} \Delta {\xi _\mathcal{A}} + \frac{1}{2}\sum\limits_{\mathcal{A},\mathcal{B} = 1}^m {{{\left. {\frac{{{\partial ^2}\mathcal{C}(\xi ,\xi \prime )}}{{\partial {{\xi '}_\mathcal{A}}\partial {{\xi '}_\mathcal{B}}}}} \right|}_{\xi \prime= \xi }}} \Delta {\xi _\mathcal{A}}\Delta {\xi _\mathcal{B}}.$$

(92)

As ${\mathcal C}$ attains its maximum value when ξ − ξ^′ = 0, we have

$${\left. {{{\partial {\mathcal C}(\xi ,{\xi \prime})} \over {\partial \xi _A\prime}}} \right|_{{\xi \prime} = \xi}} = 0,\quad{\mathcal A} = 1, \ldots ,m.$$

(93)

Let us introduce the symmetric matrix G with components

$${G_{\mathcal{A}\mathcal{B}}}(\xi ): =- {\left. {\frac{1}{{2\mathcal{C}(\xi ,\xi )}}\frac{{{\partial ^2}\mathcal{C}(\xi ,\xi \prime )}}{{\partial {{\xi '}_\mathcal{A}}\partial {{\xi '}_\mathcal{B}}}}} \right|_{\xi \prime= \xi }},\quad \mathcal{A},\mathcal{B} = 1, \ldots ,m.$$

(94)

then, the inequality (91) for the elementary cell can approximately be written as

$$\sum\limits_{\mathcal{A},\mathcal{B} = 1}^m {{G_{\mathcal{A}\mathcal{B}}}(\xi)\Delta {\xi _\mathcal{A}}\Delta {\xi _\mathcal{B}} \leq{1 \over 2}.}$$

(95)

It is interesting to find a relation between the matrix G and the Fisher matrix. One can show (see [78], Appendix B) that the matrix G is precisely equal to the reduced Fisher matrix ${\tilde \Gamma}$ given by Eq. (81).

If the components of the matrix G are constant (i.e., they are independent of the values of the intrinsic parameters ξ of the signal), the above equation defines a hyperellipsoid in m-dimensional (m is the number of the intrinsic parameters) Euclidean space ℝ^m. The m-dimensional Euclidean volume V_cell of the elementary cell given by Eq. (95) equals

$${V_{{\rm{cell}}}} = {{{{(\pi/2)}^{m/2}}} \over {\Gamma (m/2 + 1)\sqrt {\det \,{\rm{G}}}}},$$

(96)

where Γ denotes the Gamma function. We estimate the number N_cells of elementary cells by dividing the total Euclidean volume V of the m-dimensional intrinsic parameter space by the volume V_cell of one elementary cell, i.e., we have

$${N_{{\rm{cells}}}} = {V \over {{V_{{\rm{cell}}}}}}.$$

(97)

The components of the matrix G are constant for the signal h(t; A₀, ϕ₀, ξ) = A₀ cos (ϕ(t; ξ) − ϕ₀), provided the phase ϕ(t; ξ) is a linear function of the intrinsic parameters ξ.

To estimate the number of cells in the case when the components of the matrix G are not constant, i.e., when they depend on the values of the intrinsic parameters ξ, one replaces Eq. (97) by

$${N_{{\rm{cells}}}} = {{\Gamma (m/2 + 1)} \over {{{(\pi/2)}^{m/2}}}}\int\nolimits_V {\sqrt {\det {\rm{G(}}\xi)} {\rm{d}}V.}$$

(98)

This formula can be thought of as interpreting the matrix G as the metric on the parameter space. This interpretation appeared for the first time in the context of gravitational-wave data analysis in the work by Owen [102], where an analogous integral formula was proposed for the number of templates needed to perform a search for gravitational-wave signals from coalescing binaries.

The concept of number of cells was introduced in [68] and it is a generalization of the idea of an effective number of samples introduced in [46] for the case of a coalescing binary signal.

We approximate the pdf of the ${\mathcal F}$-statistic in each cell by the pdf ${p_0}({\mathcal F})$ of the ${\mathcal F}$-statistic when the parameters are known [it is given by Eq. (84)]. The values of the ${\mathcal F}$-statistic in each cell can be considered as independent random variables. The probability that ${\mathcal F}$ does not exceed the threshold ${{\mathcal F}_0}$ in a given cell is $1 - {P_F}({{\mathcal F}_0})$, where ${P_F}({{\mathcal F}_0})$ is given by Eq. (86). Consequently the probability that ${\mathcal F}$ does not exceed the threshold ${{\mathcal F}_0}$ in all the N_cells cells is ${[1 - {P_F}({{\mathcal F}_0})]^{{N_{{\rm{cells}}}}}}$. Thus, the probability $P_F^T$ that ${\mathcal F}$ exceeds ${{\mathcal F}_0}$ in one or more cells is given by

$$P_{\rm{F}}^T({\mathcal{F}_0}) = 1 - {[1 - {P_{\rm{F}}}({\mathcal{F}_0})]^{{N_{{\rm{cells}}}}}}.$$

(99)

By definition, this is the false alarm probability when the phase parameters are unknown. The number of false alarms N_F is given by

$${N_{\rm{F}}} = {N_{{\rm{cells}}}}P_{\rm{F}}^T({\mathcal{F}_0}).$$

(100)

A different approach to the calculation of the number of false alarms using the Euler characteristic of level crossings of a random field is described in [65].

It was shown (see [39]) that for any finite ${{\mathcal F}_0}$ and N_cells, Eq. (99) provides an upper bound for the false alarm probability. Also in [39] a tighter upper bound for the false alarm probability was derived by modifying a formula obtained by Mohanty [92]. The formula amounts essentially to introducing a suitable coefficient multiplying the number N_cells of cells.

4.3.3 Detection probability for unknown intrinsic parameters

When the signal is present in the data a precise calculation of the pdf of the ${\mathcal F}$-statistic is very difficult because the presence of the signal makes the data’s random process non-stationary. As a first approximation we can estimate the probability of detection of the signal when the intrinsic parameters are unknown by the probability of detection when these parameters are known [it is given by Eq. (87)]. This approximation assumes that when the signal is present the true values of the intrinsic parameters fall within the cell where the ${\mathcal F}$-statistic has a maximum. This approximation will be the better the higher the signal-to-noise ratio ρ is.

4.4 Number of templates

To search for gravitational-wave signals we evaluate the ${\mathcal F}$-statistic on a grid in parameter space. The grid has to be sufficiently fine such that the loss of signals is minimized. In order to estimate the number of points of the grid, or in other words the number of templates that we need to search for a signal, the natural quantity to study is the expectation value of the ${\mathcal F}$-statistic when the signal is present.

Thus, we assume that the data x contains the gravitational-wave signal h(t; θ) defined in Eq. (32), so x(t; θ) = h(t; θ) + n(t). The parameters θ = (a, ξ) of the signal consist of extrinsic parameters a and intrinsic parameters ξ. The data x will be correlated with the filters h_i(t; ξ^′) (i = 1,…, n) parameterized by the values ξ^′ of the intrinsic parameters. The ${\mathcal F}$-statistic can thus be written in the form [see Eq. (62)]

$$\mathcal{F}[x(t;{\rm{a,}}\xi {\rm{);}}{\xi {\prime}}] = {1 \over 2}{\rm{N}}{[x(t;{\rm{a}},\xi);{\xi {\prime}}]^{\rm{T}}} \cdot {\rm{M(}}{\xi {\prime}})^{- 1}\cdot{\rm{N}}[x(t;{\rm{a}},\xi);{\xi {\prime}}],$$

(101)

where the matrices M and N are defined in Eqs. (60). The expectation value of the ${\mathcal F}$-statistic (101) is

$${\rm{E}}[\mathcal{F}[x(t;{\rm{a}},\xi);{\xi \prime}]] = {1 \over 2}(n + {{\rm{a}}^{\rm{T}}} \cdot {\rm{Q}}(\xi, {\xi {\prime}}) \cdot {\rm{M}}{({\xi {\prime}})^{- 1}} \cdot {\rm{Q}}{(\xi, {\xi {\prime}})^{\rm{T}}} \cdot {\rm{a}}),$$

(102)

where the matrix Q is defined in Eq. (90). Let us rewrite the expectation value (102) in the following form,

$${\rm{E}}[{\mathcal F}[x(t;{\rm{a}},\xi);{\xi \prime}]] = {1 \over 2}(n + \rho {({\rm{a}},\xi)^2}{\mathcal{C}_{\rm{n}}}({\rm{a}},\xi, {\xi {\prime}})),$$

(103)

where ρ is the signal-to-noise ratio and where we have introduced the normalized correlation function ${{\mathcal C}_{\rm{n}}}$,

$${{\mathcal C}_{\rm{n}}}({\rm{a}},\xi, {\xi \prime}): = {{{{\rm{a}}^{\rm{T}}} \cdot {\rm{Q}}(\xi, {\xi \prime})\cdot{\rm{M}}{{({\xi \prime})}^{- 1}}\cdot{\rm{Q}}{{(\xi, {\xi \prime})}^{\rm{T}}}\cdot{\rm{a}}} \over {{{\rm{a}}^{\rm{T}}}\cdot{\rm{M}}(\xi)\cdot{\rm{a}}}}.$$

(104)

From the Rayleigh principle [90] it follows that the minimum value of the normalized correlation function is equal to the smallest eigenvalue of the matrix M(ξ)⁻¹ · Q(ξ, ξ^′) · M(ξ^′)⁻¹ · Q(ξ, ξ^′)^T, whereas the maximum value is given by its largest eigenvalue. We define the reduced correlation function ${\mathcal C}$ as

$$\mathrm{C}(\xi, {\xi {\prime}}): = {1 \over 2}{\rm{Tr}}({\rm{M}}{(\xi)^{- 1}} \cdot {\rm{Q}}(\xi, {\xi {\prime}}) \cdot {\rm{M}}{({\xi {\prime}})^{- 1}} \cdot {\rm{Q}}{(\xi, {\xi {\prime}})^{\rm{T}}}).$$

(105)

As the trace of a matrix equals the sum of its eigenvalues, the reduced correlation function ${\mathcal C}$ is equal to the average of the eigenvalues of the normalized correlation function ${{\mathcal C}_{\rm{n}}}$. In this sense we can think of the reduced correlation function as an “average” of the normalized correlation function. The advantage of the reduced correlation function is that it depends only on the intrinsic parameters ξ, and thus is suitable for studying the number of grid points on which the ${\mathcal F}$-statistic needs to be evaluated. We also note that the normalized correlation function ${\mathcal C}$ precisely coincides with the autocovariance function ${\mathcal C}$ of the ${\mathcal F}$-statistic given by Eq. (89).

As in the calculation of the number of cells in order to estimate the number of templates we perform a Taylor expansion of ${\mathcal C}$ up to second order terms around the true values of the parameters, and we obtain an equation analogous to Eq. (95),

$$\sum\limits_{\mathcal{A},\mathcal{B} = 1}^m {{G_{\mathcal{A}\mathcal{B}}}\Delta {\xi _\mathcal{A}}\Delta {\xi _\mathcal{B}} = 1 - {C_0},}$$

(106)

where G is given by Eq. (94). By arguments identical to those in deriving the formula for the number of cells we arrive at the following formula for the number of templates:

$${N_t} = {1 \over {{{(1 - {C_0})}^{m/2}}}}{{\Gamma (m/2 + 1)} \over {{\pi ^{m/2}}}}\int\nolimits_V {\sqrt {\det G(\xi)} {\rm d}V.}$$

(107)

When C₀ = 1/2 the above formula coincides with the formula for the number N_cells of cells, Eq. (98). Here we would like to place the templates sufficiently closely so that the loss of signals is minimized. Thus 1 − C₀ needs to be chosen sufficiently small. The formula (107) for the number of templates assumes that the templates are placed in the centers of hyperspheres and that the hyperspheres fill the parameter space without holes. In order to have a tiling of the parameter space without holes we can place the templates in the centers of hypercubes, which are inscribed in the hyperspheres. Then the formula for the number of templates reads

$${N_t} = {1 \over {{{(1 - {C_0})}^{m/2}}}}{{{m^{m/2}}} \over {{2^m}}}\int\nolimits_V {\sqrt {\det {\rm{G}}(\xi)} {\rm{d}}V.}$$

(108)

For the case of the signal given by Eq. (34) our formula for the number of templates is equivalent to the original formula derived by Owen [102]. Owen [102] has also introduced a geometric approach to the problem of template placement involving the identification of the Fisher matrix with a metric on the parameter space. An early study of the template placement for the case of coalescing binaries can be found in [121, 45, 26]. Applications of the geometric approach of Owen to the case of spinning neutron stars and supernova bursts are given in [33, 16].

4.4.1 Covering problem

The problem of how to cover the parameter space with the smallest possible number of templates, such that no point in the parameter space lies further away from a grid point than a certain distance, is known in mathematical literature as the covering problem [38]. This was first studied in the context of gravitational-wave data analysis by Prix [ ]. The maximum distance of any point to the next grid point is called the covering radius R. An important class of coverings are lattice coverings. We define a lattice in m-dimensional Euclidean space ℝ^m to be the set of points including 0 such that if u and v are lattice points, then also u + v and u − v are lattice points. The basic building block of a lattice is called the fundamental region. A lattice covering is a covering of ℝ^m by spheres of covering radius R, where the centers of the spheres form a lattice. The most important quantity of a covering is its thickness Θ defined as

$$\Theta: = {{{\rm{volume}}\,{\rm{of}}\,{\rm{one}}\,m-{\rm{dimensional}}\,{\rm{sphere}}} \over {{\rm{volume}}\,{\rm{of}}\,{\rm{the}}\,{\rm{fundamental}}\,{\rm{region}}}}.$$

(109)

In the case of a two-dimensional Euclidean space the best covering is the hexagonal covering and its thickness ≃ 1.21. For dimensions higher than 2 the best covering is not known. However, we know the best lattice covering for dimensions m ≤ 23. These are A*_m lattices, which have thicknesses ${\Theta _{A_m^{\ast}}}$ equal to

$${\Theta _{A_m^{\ast}}} = {V_m}\sqrt {m + 1} {\left({{{m(m + 2)} \over {12(m + 1)}}} \right)^{m/2}},$$

(110)

where V_m is the volume of the m-dimensional sphere of unit radius. The advantage of an A*_m lattice over the hypercubic lattice grows exponentially with the number of dimensions.

For the case of gravitational-wave signals from spinning neutron stars a 3-dimensional grid was constructed [18]. It consists of prisms with hexagonal bases. Its thickness is around 1.84, which is much better than the cubic grid with a thickness of approximately 2.72. It is worse than the best 3-dimensional lattice covering, which has a thickness of around 1.46.

In [19] a grid was constructed in the 4-dimensional parameter space spanned by frequency, frequency derivative, and sky position of the source, for the case of an almost monochromatic gravitational-wave signal originating from a spinning neutron star. The starting point of the construction was an A*₄ lattice of thickness ≃ 1.77. The grid was then constrained so that the nodes of the grid coincide with Fourier frequencies. This allowed the use of a fast Fourier transform (FFT) to evaluate the maximum-likelihood ${\mathcal F}$-statistic efficiently (see Section 4.6.2). The resulting lattice is only 20% thicker than the optimal A*₄ lattice.

Efficient 2-dimensional banks of templates suitable for directed searches (in which one assumes that the position of the gravitational-wave source in the sky is known, but one does not assume that the wave’s frequency and its derivative are a priori known) were constructed in [104]. All grids found in [104] enable usage of the FFT algorithm in the computation of the ${\mathcal F}$-statistic; they have thicknesses 0.1–16% larger than the thickness of the optimal 2-dimensional hexagonal covering. In the construction of grids the dependence on the choice of the position of the observational interval with respect to the origin of time axis was employed. Also the usage of the FFT algorithms with nonstandard frequency resolutions achieved by zero padding or folding the data was discussed.

The above template placement constructions are based on a Fisher matrix with constant coefficients, i.e., they assume that the parameter manifold is flat. The generalization to curved Riemannian parameter manifolds is difficult. An interesting idea to overcome this problem is to use stochastic template banks where a grid in the parameter space is randomly generated by some algorithm [89, 57, 86, 119].

4.5 Suboptimal filtering

To extract gravitational-wave signals from the detector’s noise one very often uses filters that are not optimal. We may have to choose an approximate, suboptimal filter because we do not know the exact form of the signal (this is almost always the case in practice) or in order to reduce the computational cost and to simplify the analysis. In the case of the signal of the form given in Eq. (32) the most natural and simplest way to proceed is to use as detection statistic the ${\mathcal F}$ statistic where the filters h′_k(t; ζ) (k = 1, …, n) are the approximate ones instead of the optimal ones h_k(t; ξ) (k = 1,…, n) matched to the signal. In general the functions h′_k(t; ζ) will be different from the functions h_k(t; ξ) used in optimal filtering, and also the set of parameters ζ will be different from the set of parameters ξ in optimal filters. We call this procedure the suboptimal filtering and we denote the suboptimal statistic by ${{\mathcal F}_{\rm{s}}}$. It is defined as [see Eq. (62)]

$${\mathcal{F}_{\rm{s}}}[x;\zeta ]: = {1 \over 2}{{\rm{N}}_{\rm{s}}}{[x;\zeta ]^{\rm T}} \cdot {{\rm{M}}_{\rm{s}}}{(\zeta)^{- 1}} \cdot {{\rm{N}}_{\rm{s}}}[x;\zeta ],$$

(111)

where the data-dependent n × 1 column matrix N_s and the square n × n matrix M_s have components [see Eq. (60)]

$${N_{{\text{s}}\;i}}[x;\zeta ]: = (x(t)|{h'_i}(t;\zeta )),\quad {M_{{\text{s}}ij}}(\zeta ): = ({h'_i}(t;\zeta )|{h'_j}(t;\zeta )),\quad i,j = 1, \ldots ,n.$$

(112)

We need a measure of how well a given suboptimal filter performs. To find such a measure we calculate the expectation value of the suboptimal statistic ${{\mathcal F}_{\rm{s}}}$ in the case where the data contains the gravitational-wave signal, i.e., when x(t; a, ξ) = n(t) + h(t; a, ξ). We get

$${\rm{E}}[{{\mathcal{F}}_s}[x(t;{\rm{a}},\xi);\zeta ]] = {1 \over 2}(n + {{\rm{a}}^{\rm T}} \cdot {\rm{Q}} _{\rm{s}}(\xi, \zeta) \cdot {{\rm{M}}_{\rm{s}}}{(\zeta)^{- 1}} \cdot {\rm Q}_{\rm{s}}{(\xi, \zeta)^{\rm{T}}} \cdot \rm{a}),$$

(113)

where we have introduced the matrix Q_s with components

$${Q_{{\text{s}}ij}}(\xi ,\zeta ): = ({h_i}(t;\xi )|{h'_j}(t;\zeta )),\quad i,j = 1, \ldots ,n.$$

(114)

Let us rewrite the expectation value (113) in the following form,

$${\rm{E}}[{{\mathcal F}_{\rm{s}}}[x(t;{\rm{a}},\xi);\zeta ]] = {1 \over 2}\left({n + \rho {{({\rm{a}},\xi)}^2}{{{{\rm{a}}^T} \cdot {{\rm{Q}}_{\rm{s}}}(\xi, \zeta) \cdot {{\rm{M}}_{\rm{s}}}{{(\zeta)}^{- 1}} \cdot {{\rm{Q}}_{\rm{s}}}{{(\xi, \zeta)}^T} \cdot {\rm{a}}} \over {{{\rm{a}}^{\rm{T}}} \cdot M(\xi) \cdot {\rm{a}}}}} \right),$$

(115)

where ρ is the optimal signal-to-noise ratio [given in Eq. (64)]. This expectation value reaches its maximum equal to (n + ρ²)/2 when the filter is perfectly matched to the signal. Therefore, a natural measure of the performance of a suboptimal filter is the quantity FF defined by

$$\text{FF}(\xi ): = \mathop {\max }\limits_{(a,\varsigma )} \sqrt {\frac{{{\text{a}^\text{T}} \cdot {\text{Q}_\text{s}}(\xi ,\varsigma ) \cdot {\text{M}_\text{s}}{{(\varsigma )}^{ - 1}} \cdot {\text{Q}_\text{s}}{{(\xi ,\varsigma )}^\text{T}} \cdot \text{a}}}{{{\text{a}^\text{T}} \cdot \text{M}(\xi ) \cdot \text{a}}}} .$$

(116)

We call the quantity FF the generalized fitting factor. From the Rayleigh principle, it follows that the generalized fitting factor is the maximum of the largest eigenvalue of the matrix M(ξ)⁻¹ · Q_s(ξ, ζ) · M_s(ζ)⁻¹ · Q_s(ξ, ζ)^T over the intrinsic parameters of the signal.

In the case of a gravitational-wave signal given by

$$s(t;{A_0},\xi) = {A_0}h(t;\xi),$$

(117)

the generalized fitting factor defined above reduces to the fitting factor introduced by Apostolatos [13]:

$$\mathrm{FF}(\xi) = \max\limits_\zeta {{(h(t;\xi)\vert {h{\prime}}(t;\zeta))} \over {\sqrt {(h(t;\xi)\vert h(t;\xi))} \sqrt {({h{\prime}}(t;\zeta)\vert {h{\prime}}(t;\zeta))}}}.$$

(118)

The fitting factor is the ratio of the maximal signal-to-noise ratio that can be achieved with suboptimal filtering to the signal-to-noise ratio obtained when we use a perfectly matched, optimal filter. We note that for the signal given by Eq. (117), FF is independent of the value of the amplitude A₀.

For the case of a signal of the form

$$s(t;{A_0},{\phi _0},\xi) = {A_0}\cos (\phi (t;\xi) + {\phi _0}),$$

(119)

where ϕ₀ is a constant phase, the maximum over ϕ₀ in Eq. (118) can be obtained analytically. Moreover, assuming that over the bandwidth of the signal the spectral density of the noise is constant and that over the observation time cosϕ(t; ξ) oscillates rapidly, the fitting factor is approximately given by

$${\rm{FF}}(\xi) \cong \max\limits_\zeta {\left[ {{{\left({\int\nolimits_0^{{T_0}} {\cos (\phi (t;\xi) - {\phi {\prime}}(t;\zeta))dt}} \right)}^2} + {{\left({\int\nolimits_0^{{T_0}} {\sin (\phi (t;\xi) - {\phi {\prime}}(t;\zeta))dt}} \right)}^2}} \right]^{1/2}}.$$

(120)

In designing suboptimal filters one faces the issue of how small a fitting factor one can accept. A popular rule of thumb is accepting FF = 0. 97. Assuming that the amplitude of the signal and consequently the signal-to-noise ratio decreases inversely proportionally to the distance from the source this corresponds to 10% loss of the signals that would be detected by a matched filter.

Proposals for good suboptimal (search) templates for the case of coalescing binaries are given in [35, 134] and for the case-spinning neutron stars in [65, 18].

4.6 Algorithms to calculate the ${\mathcal F}$-statistic

4.6.1 The two-step procedure

In order to detect signals we search for threshold crossings of the ${\mathcal F}$-statistic over the intrinsic parameter space. Once we have a threshold crossing we need to find the precise location of the maximum of ${\mathcal F}$ in order to estimate accurately the parameters of the signal. A satisfactory procedure is the two-step procedure. The first step is a coarse search where we evaluate ${\mathcal F}$ on a coarse grid in parameter space and locate threshold crossings. The second step, called a fine search, is a refinement around the region of parameter space where the maximum identified by the coarse search is located.

There are two methods to perform the fine search. One is to refine the grid around the threshold crossing found by the coarse search [94, 92, 134, 127], and the other is to use an optimization routine to find the maximum of ${\mathcal F}$ [65, 78]. As initial values to the optimization routine we input the values of the parameters found by the coarse search. There are many maximization algorithms available. One useful method is the Nelder-Mead algorithm [79], which does not require computation of the derivatives of the function being maximized.

4.6.2 Evaluation of the ${\mathcal F}$-statistic

Usually the grid in parameter space is very large and it is important to calculate the optimum statistic as efficiently as possible. In special cases the ${\mathcal F}$-statistic given by Eq. (62) can be further simplified. For example, in the case of coalescing binaries ${\mathcal F}$ can be expressed in terms of convolutions that depend on the difference between the time-of-arrival (TOA) of the signal and the TOA parameter of the filter. Such convolutions can be efficiently computed using FFTs. For continuous sources, like gravitational waves from rotating neutron stars observed by ground-based detectors [65] or gravitational waves form stellar mass binaries observed by space-borne detectors [78], the detection statistic ${\mathcal F}$ involves integrals of the general form

$$\int\nolimits_0^{{T_0}} {x(t)m(t;\omega, \tilde \xi)} \exp (i\omega {\phi _{\bmod}}(t;\tilde \xi))\exp (i\omega t){\rm{d}}t,$$

(121)

where ${\tilde \xi}$ are the intrinsic parameters excluding the frequency parameter ω, m is the amplitude modulation function, and ωϕ_mod the phase modulation function. The amplitude modulation function is slowly varying compared to the exponential terms in the integral (121). We see that the integral (121) can be interpreted as a Fourier transform (and computed efficiently with an FFT), if ϕ_mod = 0 and if m does not depend on the frequency ω. In the long-wavelength approximation the amplitude function m does not depend on the frequency. In this case, Eq. (121) can be converted to a Fourier transform by introducing a new time variable t_b [124],

$${t_{\rm{b}}}(t;\tilde \xi): = t + {\phi _{\bmod}}(t;\tilde \xi).$$

(122)

Thus, in order to compute the integral (121), for each set of the intrinsic parameters ${\tilde \xi}$ we multiply the data by the amplitude modulation function m, resample according to Eq. (122), and perform the FFT. In the case of LISA detector data when the amplitude modulation m depends on frequency we can divide the data into several band-passed data sets, choosing the bandwidth for each set to be sufficiently small so that the change of m exp(iωϕ_mod) is small over the band. In the integral (121) we can then use as the value of the frequency in the amplitude and phase modulation function the maximum frequency of the band of the signal (see [78] for details).

4.7 Accuracy of parameter estimation

4.7.1 Fisher-matrix-based assessments

Fisher matrix has been extensively used to assess the accuracy of estimation of astrophysically-interesting parameters of different gravitational-wave signals. For ground-based interferometric detectors, the first calculations of the Fisher matrix concerned gravitational-wave signals from inspiralling compact binaries (made of neutron stars or black holes) in the leading-order quadrupole approximation [50, 76, 63] and from quasi-normal modes of Kerr black hole [48].

Cutler and Flanagan [41] initiated the study of the implications of the higher-order post-Newtonian (PN) phasing formula as applied to the parameter estimation of inspiralling binary signals. They used the 1.5PN phasing formula to investigate the problem of parameter estimation, both for spinning and non-spinning binaries, and examined the effect of the spin-orbit coupling on the estimation of parameters. The effect of the 2PN phasing formula was analyzed independently by Poisson and Will [106] and Królak, Kokkotas and Schäfer [75]. In both cases the focus was to understand the leading-order spin-spin coupling term appearing at the 2PN level when the spins were aligned perpendicularly to the orbital plane. Compared to [75], [106] also included a priori information about the magnitude of the spin parameters, which then leads to a reduction in the rms errors in the estimation of mass parameters. The case of a 3.5PN phasing formula was studied in detail by Arun et al. [17]. Inclusion of 3.5PN effects leads to an improved estimate of the binary parameters. Improvements are relatively smaller for lighter binaries. More recently the Fisher matrix was employed to assess the errors in estimating the parameters of nonspinning black-hole binaries using the complete inspiral-merger-ring-down waveforms [7].

Various authors have investigated the accuracy with which the LISA detector can determine binary parameters including spin effects. Cutler [40] determined LISA’s angular resolution and evaluated the errors of the binary masses and distance considering spins aligned or anti-aligned with the orbital angular momentum. Hughes [60] investigated the accuracy with which the redshift can be estimated (if the cosmological parameters are derived independently), and considered the black-hole ring-down phase in addition to the inspiralling signal. Seto [128] included the effect of finite armlength (going beyond the long wavelength approximation) and found that the accuracy of the distance determination and angular resolution improve. This happens because the response of the instrument when the armlength is finite depends strongly on the location of the source, which is tightly correlated with the distance and the direction of the orbital angular momentum. Vec-chio [140] provided the first estimate of parameters for precessing binaries when only one of the two supermassive black holes carries spin. He showed that modulational effects decorrelate the binary parameters to some extent, resulting in a better estimation of the parameters compared to the case when spins are aligned or antialigned with orbital angular momentum. Hughes and Menou [61] studied a class of binaries, which they called “golden binaries,” for which the inspiral and ring-down phases could be observed with good enough precision to carry out valuable tests of strong-field gravity. Berti, Buonanno and Will [29] have shown that inclusion of non-precessing spin-orbit and spin-spin terms in the gravitational-wave phasing generally reduces the accuracy with which the parameters of the binary can be estimated. This is not surprising, since the parameters are highly correlated, and adding parameters effectively dilutes the available information.

Extensive study of accuracy of parameter estimation for continuous gravitational-wave signals from spinning neutron stars was performed in [64]. In [129] Seto used the Fisher matrix to study the possibility of determining distances to rapidly rotating isolated neutron stars by measuring the curvature of the wave fronts.

4.7.2 Comparison with the Cramèr-Rao bound

In order to test the performance of the maximization method of the ${\mathcal F}$-statistic it is useful to perform Monte Carlo simulations of the parameter estimation and compare the simulated variances of the estimators with the variances calculated from the Fisher matrix. Such simulations were performed for various gravitational-wave signals [73, 26, 65, 36]. In these simulations one observes that, above a certain signal-to-noise ratio, called the threshold signal-to-noise ratio, the results of the Monte Carlo simulations agree very well with the calculations of the rms errors from the inverse of the Fisher matrix. However, below the threshold signal-to-noise ratio they differ by a large factor. This threshold effect is well known in signal processing [139]. There exist more refined theoretical bounds on the rms errors that explain this effect, and they were studied in the context of the gravitational-wave signals from coalescing binaries [98].

Use of the Fisher matrix in the assessment of accuracy of the parameter estimation has been critically examined in [138], where a criterion has been established for the signal-to-noise ratio above which the inverse of the Fisher matrix approximates well covariances of the parameter estimators. In [148, 142] the errors of ML estimators of parameters of gravitational-wave signals from nonspinning black-hole binaries were calculated analytically using a power expansion of the bias and the covariance matrix in inverse powers of the signal-to-noise ratio. The first-order term in this covariance matrix expansion is the inverse of the Fisher information matrix. The use of higher-order derivatives of the likelihood function in these expansions makes the errors prediction sensitive to the secondary lobes of the pdf of the ML estimators. Conditions for the validity of the Cramèr-Rao lower bound are discussed in [142] as well, and some new features in regions of the parameter space so far not explored are predicted (e.g., that the bias can become the most important contributor to the parameters errors for high-mass systems with masses 200 M_⊙ and above).

There exists a simple model that explains the deviations from the covariance matrix and reproduces well the results of the Monte Carlo simulations (see also [25]). The model makes use of the concept of the elementary cell of the parameter space that we introduced in Section 4.3.2. The calculation given below is a generalization of the calculation of the rms error for the case of a monochromatic signal given by Rife and Boorstyn [116].

When the values of parameters of the template that correspond to the maximum of the functional ${\mathcal F}$ fall within the cell in the parameter space where the signal is present, the rms error is satisfactorily approximated by the inverse of the Fisher matrix. However, sometimes, as a result of noise, the global maximum is in the cell where there is no signal. We then say that an outlier has occurred. In the simplest case we can assume that the probability density of the values of the outliers is uniform over the search interval of a parameter, and then the rms error is given by

$$\sigma _{{\rm{out}}}^2 = {{{\Delta ^2}} \over {12}},$$

(123)

where Δ is the length of the search interval for a given parameter. The probability that an outlier occurs will be higher the lower the signal-to-noise ratio is. Let q be the probability that an outlier occurs. Then the total variance σ² f the estimator of a parameter is the weighted sum of the two errors

$${\sigma ^2} = \sigma _{{\rm{out}}}^2q + \sigma _{{\rm{CR}}}^2(1 - q),$$

(124)

where σ_CR is the rms errors calculated from the covariance matrix for a given parameter. One can show [65] that the probability q can be approximated by the following formula:

$$q = 1 - \int\nolimits_0^\infty {{p_1}(\rho, \mathcal{F}){{\left({\int\nolimits_0^\mathcal{F} {{p_0}(y){\rm{d}}y}} \right)}^{{N_{{\rm{cells}}}} - 1}}{\rm{d}}\mathcal{F},}$$

(125)

where p₀ and p₁ are the pdfs of the ${\mathcal F}$-statistic (for known intrinsic parameters) when the signal is absent or present in data, respectively [they are given by Eqs. (84) and (85)], and where N_cells is the number of cells in the intrinsic parameter space. Eq. (125) is in good but not perfect agreement with the rms errors obtained from the Monte Carlo simulations (see [65]). There are clearly other reasons for deviations from the Cramèr-Rao bound as well. One important effect (see [98]) is that the functional ${\mathcal F}$ has many local subsidiary maxima close to the global one. Thus, for a low signal-to-noise ratio the noise may promote the subsidiary maximum to a global one.

4.8 Upper limits

Detection of a signal is signified by a large value of the ${\mathcal F}$-statistic that is unlikely to arise from the noise-only distribution. If instead the value of ${\mathcal F}$ is consistent with pure noise with high probability we can place an upper limit on the strength of the signal. One way of doing this is to take the loudest event obtained in the search and solve the equation

$${P_{\rm{D}}}({\rho _{{\rm{UL}}}},{{\mathcal F}_{\rm{L}}}) = \beta$$

(126)

for signal-to-noise ratio ρ_UL, where P_D is the detection probability given by Eq. (87), ${{\mathcal F}_{\rm{L}}}$ is the value of the ${\mathcal F}$-statistic corresponding to the loudest event, and β is a chosen confidence [23, 2]. Then ρ_UL is the desired upper limit with confidence β.

When gravitational-wave data do not conform to a Gaussian probability density assumed in Eq. (87), a more accurate upper limit can be obtained by injecting the signals into the detector’s data and thereby estimating the probability of detection P_D [4].

5 Network of Detectors

Several gravitational-wave detectors can observe gravitational waves from the same source. For example a network of bar detectors can observe a gravitational-wave burst from the same supernova explosion, or a network of laser interferometers can detect the inspiral of the same compact binary system. The space-borne LISA detector can be considered as a network of three detectors that can make three independent measurements of the same gravitational-wave signal. Simultaneous observations are also possible among different types of detectors. For example, a search for supernova bursts can be performed simultaneously by resonant and interferometric detectors [21].

Let us consider a network consisting of N gravitational-wave detectors and let us denote by x_I the data collected by the Ith detector (I = 1, …, N). We assume that noises in all detectors are additive, so the data x_I is a sum of the noise n_I in the Ith detector and eventually a gravitational-wave signal h_I registered by the Ith detector,

$${x_I}(t) = {n_I}(t) + {h_I}(t),\,I = 1, \ldots, N.$$

(127)

It is convenient to collect all the data streams, all the noises, and all the gravitational-wave signals into column N × 1 matrices denoted respectively by x, n, and h,

$${\rm{x}}(t): = \left({\begin{array}{*{20}c} {{x_1}(t)} \\ \vdots \\ {{x_N}(t)} \\ \end{array}} \right),{\rm{n}}(t): = \left({\begin{array}{*{20}c} {{n_1}(t)} \\ \vdots \\ {{n_N}(t)} \\ \end{array}} \right),\,{\rm{h}}(t): = \left({\begin{array}{*{20}c} {{h_1}(t)} \\ \vdots \\ {{h_N}(t)} \\ \end{array}} \right),$$

(128)

then Eqs. (127) can shortly be written as

$${\rm{x}}(t) = {\rm{n}}(t) + {\rm{h}}(t).$$

(129)

if additionally all detectors’ noises are stationary, Gaussian, and continuous-in-time random processes with zero means, the network log likelihood function is given by

$$\log \Lambda [{\rm{x}}] = ({\rm{x}} + {\rm{h}}) - {1 \over 2}({\rm{h}} + {\rm{h}}),$$

(130)

where the scalar product (· ∣ ·) is defined by

$$({\rm{x}}\vert {\rm{y}}): = 4\Re \int\nolimits_0^\infty {{\rm{\tilde x}}{{(f)}^{\rm{T}}} \cdot {{\rm{S}}_n}{{(f)}^{- 1}} \cdot {\rm{\tilde y}}} {(f)^{\ast}}{\rm{d}}f.$$

(131)

Here S_n is the one-sided cross spectral density matrix of the noises of the detector network, which is defined by (here E denotes the expectation value)

$${\rm{E}}[{\rm{\tilde n}}(f) \cdot {{\rm{\tilde n}}^{\ast}}{({f{\prime}})^{\rm{T}}}] = {1 \over 2}\delta (f - {f{\prime}}){{\rm{S}}_n}(\vert f\vert).$$

(132)

The analysis is greatly simplified if the cross spectrum matrix S_n is diagonal. This means that the noises in various detectors are uncorrelated. This is the case when the detectors of the network are in widely separated locations, like, for example, the two LIGO detectors. However, this assumption is not always satisfied. An important case is the LISA detector where the noises of the three independent responses are correlated. Nevertheless for the case of LISA one can find a set of three combinations for which the noises are uncorrelated [108, 112]. When the cross spectrum matrix is diagonal the network log likelihood function is just the sum of the log likelihood functions for each detector.

Derivation of the likelihood function for an arbitrary network of detectors can be found in [49]. Applications of optimal filtering for observations of gravitational-wave signals from coalescing binaries by networks of ground-based detectors are given in [63, 41, 62, 103], and for the case of stellar-mass binaries observed by LISA space-borne detector in [78, 118]. The single-detector ${\mathcal F}$-statistic for nearly monochromatic gravitational waves from spinning neutron stars was generalized to the case of a network of detectors (also with time-varying noise curves) in [42] (in this work the ${\mathcal F}$-statistic was also generalized from the usual single-source case to the case of a collection of known sources). The reduced Fisher matrix [defined in Eq. (81)] for the case of a network of interferometers observing spinning neutron stars has been derived and studied in [110].

Network searches for gravitational-wave burst signals of unknown shape are often based on maximization of the network likelihood function over each sample of the unknown polarization waveforms h₊ and h_× and over sky positions of the source [52, 95]. A least-squares-fit solution for the estimation of the sky location of the source and the polarization waveforms by a network of three detectors for the case of a broadband burst was obtained in [56].

There is also another important method for analyzing the data from a network of detectors — the search for coincidences of events among detectors. This analysis is particularly important when we search for supernova bursts, the waveforms of which are not very well known. Such signals can be easily mimicked by the non-Gaussian behavior of the detector noise. The idea is to filter the data optimally in each of the detectors and obtain candidate events. Then one compares parameters of candidate events, like, for example, times of arrivals of the bursts, among the detectors in the network. This method is widely used in the search for supernovae by networks of bar detectors [24]. A new geometric coincident algorithm of combining the data from a network of detectors was proposed in [117]. This algorithm employs the covariances between signal’s parameters in such a way that it associates with each candidate event an ellipsoidal region in parameter space defined by the covariance matrix. Events from different detectors are deemed to be in coincidence if their ellipsoids have a nonzero overlap. The coincidence and the coherent strategies of multidetector detection of gravitational-wave signals from inspiralling compact binaries have been compared in [96, 97, 32]. [96] considered detectors in pairs located in the same site and [97] pairs of detectors at geographically separated sites. The case of three detectors (like the network of two LIGO detectors and the Virgo detector) has been considered in detail in [32], where it was demonstrated that the hierarchical coherent pipeline on Gaussian data has a better performance than the pipeline with just the coincident stage.

A general framework for studying the effectiveness of networks of interferometric gravitational-wave detectors has been proposed in [125]. Using this framework it was shown that adding a fourth detector to the existing network of LIGO/VIRGO detectors can dramatically increase, by a factor of 2 to 4, the detected event rate by allowing coherent data analysis to reduce the spurious instrumental coincident background.

6 Non-stationary, Non-Gaussian, and Non-linear Data

Eqs. (61) and (62) provide maximum likelihood estimators only when the noise in which the signal is buried is Gaussian. There are general theorems in statistics indicating that the Gaussian noise is ubiquitous. One is the central limit theorem, which states that the mean of any set of variables with any distribution having a finite mean and variance tends to the normal distribution. The other comes from the information theory and says that the probability distribution of a random variable with a given mean and variance, which has the maximum entropy (minimum information) is the Gaussian distribution. Nevertheless, analysis of the data from gravitational-wave detectors shows that the noise in the detector may be non-Gaussian (see, e.g., Figure 6 in [22]). The noise in the detector may also be a non-linear and a non-stationary random process.

The maximum likelihood method does not require that the noise in the detector be Gaussian or stationary. However, in order to derive the optimum statistic and calculate the Fisher matrix we need to know the statistical properties of the data. The probability distribution of the data may be complicated, and the derivation of the optimum statistic, the calculation of the Fisher matrix components and the false alarm probabilities may be impractical. However, there is one important result that we have already mentioned. The matched-filter, which is optimal for the Gaussian case is also a linear filter that gives maximum signal-to-noise ratio no matter what the distribution of the data. Monte Carlo simulations performed by Finn [49] for the case of a network of detectors indicate that the performance of matched-filtering (i.e., the maximum likelihood method for Gaussian noise) is satisfactory for the case of non-Gaussian and stationary noise.

Allen et al. [10, 11] derived an optimal (in the Neyman-Pearson sense, for weak signals) signal processing strategy, when the detector noise is non-Gaussian and exhibits tail terms. This strategy is robust, meaning that it is close to optimal for Gaussian noise but far less sensitive than conventional methods to the excess large events that form the tail of the distribution. This strategy is based on a locally optimal test [71] that amounts to comparing a first non-zero derivative

$${\Lambda _n}[x] = {\left. {{{{d^n}\Lambda [x\vert \epsilon]} \over {d\epsilon^n}}} \right\vert _{\epsilon = 0}}$$

(133)

of the likelihood ratio with respect to the amplitude of the signal with a threshold instead of the likelihood ratio itself.

The non-stationarity in the case of Gaussian and uncorrelated noise can be easily incorporated into matched filtering (see Appendix C of [1]). Let us assume that a noise sample n_l in the data has a Gaussian pdf with a variance $\sigma _l^2$ and zero mean (l = 1,…, N, where N is the number of data points). Different noise samples may have distributions with different variances. We also assume that the noise samples are uncorrelated, then the autocorrelation function K(l, l^′) of the noise is given by [see Eq. (39)]

$$K(l,{l{\prime}}) = \sigma _l^2{\delta _{l{l{\prime}}}},$$

(134)

where δ_ll′ is the Kronecker delta function. In the case of a known signal h_l and additive noise the optimal filter q_l is the solution of the following equation [which is a discrete version of the integral Eq. (41)]:

$${h_l} = \sum\limits_{{l{\prime}} = 1}^N {K(l,{l{\prime}})ql{\prime}.}$$

(135)

Thus, we have the following equation for the filter q_l:

$${q_l} = {{{h_l}} \over {\sigma _l^2}},$$

(136)

and the following expression for the log likelihood ratio:

$$\ln \Lambda [x] = \sum\limits_{l = 1}^N {{{{h_l}{x_l}} \over {\sigma _l^2}} - {1 \over 2}} \sum\limits_{l = 1}^N {{{h_l^2} \over {\sigma _l^2}}.}$$

(137)

Thus, we see that for non-stationary, uncorrelated Gaussian noise the optimal processing is identical to matched filtering for a known signal in stationary Gaussian noise, except that we divide both the data x_l and the signal h_l by time-varying standard deviation of the noise. This may be thought of as a special case of whitening the data and then correlating it using a whitened filter.

In the remaining part of this section we review some statistical tests and methods to detect non-Gaussianity, non-stationarity, and non-linearity in the data. A classical test for a sequence of data to be Gaussian is the Kolmogorov-Smirnov test [37]. It calculates the maximum distance between the cumulative distribution of the data and that of a normal distribution, and assesses the significance of the distance. A similar test is the Lillifors test [37], but it adjusts for the fact that the parameters of the normal distribution are estimated from the data rather than specified in advance. Another test is the Jarque-Bera test [69], which determines whether sample skewness and kurtosis are unusually different from their Gaussian values.

A useful test to detect outliers in the data is Grubbs’ test [55]. This test assumes that the data has an underlying Gaussian probability distribution but it is corrupted by some disturbances. Grubbs’ test detects outliers iteratively. Outliers are removed one by one and the test is iterated until no outliers are detected. Grubbs’ test is a test of the null hypothesis:

• 1. H₀: There are no outliers in the data set x_k.

  against the alternate hypothesis:

• 2. H₁: There is at least one outlier in the data set x_k.

The Grubbs’ test statistic is the largest absolute deviation from the sample mean in units of the sample standard deviation, so it is defined as

$$G = {{{{\max}_k}\vert {x_k} - \mu \vert} \over \sigma},$$

(138)

where μ and σ denote the sample mean and the sample standard deviation, respectively. The hypothesis of no outliers is rejected if

$$G > {{n - 1} \over {\sqrt n}}\sqrt {{{t_{\alpha/2(2n),n - 2}^2} \over {n - 2 + t_{\alpha/(2n),n - 2}^2}},}$$

(139)

where t_{α/(2n),n−2} denotes the critical value of the t-distribution with n − 2 degrees of freedom and a significance level of α/(2n).

Grubbs’ test has been used to identify outliers in the search of Virgo data for gravitational-wave signals from the Vela pulsar [1]. A test to discriminate spurious events due to non-stationarity and non-Gaussianity of the data from genuine gravitational-wave signals has been developed by Allen [9]. This test, called the χ² time-frequency discriminator, is applicable to the case of broadband signals, such as those coming from compact coalescing binaries.

Let now x_k and u_l be two discrete-in-time random processes (−∞ < k, l < ∞) and let u_l be independent and identically distributed (i.i.d.) random variables. We call the process x_k linear if it can be represented by

$${x_k} = \sum\limits_{l = 0}^N {{a_l}{u_{k - l}},}$$

(140)

where a_l are constant coefficients. If u_l is Gaussian (non-Gaussian), we say that x_l is linear Gaussian (non-Gaussian). In order to test for linearity and Gaussianity we examine the third-order cumulants of the data. The third-order cumulant C_kl of a zero mean stationary process is defined by

$${C_{kl}}: = {\rm{E}}[{x_m}{x_{m + k}}{x_{m + l}}].$$

(141)

The bispectrum S₂(f₁, f₂) is the two-dimensional Fourier transform of C_kl. The bicoherence is defined as

$$B({f_1},{f_2}): = {{{S_2}({f_1},{f_2})} \over {S({f_1} + {f_2})S({f_1})S({f_2})}},$$

(142)

where S(f) is the spectral density of the process x_k. If the process is Gaussian, then its bispectrum and consequently its bicoherence is zero. One can easily show that if the process is linear then its bicoherence is constant. Thus, if the bispectrum is not zero, then the process is non-Gaussian; if the bicoherence is not constant then the process is also non-linear. Consequently we have the following hypothesis testing problems:

• 1. H₁: The bispectrum of x_k is nonzero.

• 2. H₀: The bispectrum of x_k is zero.

If Hypothesis 1 holds, we can test for linearity, that is, we have a second hypothesis testing problem:

• 3. H′₁: The bicoherence of x_k is not constant.

• 4. H″₁: The bicoherence of x_k is a constant.

If Hypothesis 4 holds, the process is linear.

Using the above tests we can detect non-Gaussianity and, if the process is non-Gaussian, non-linearity of the process. The distribution of the test statistic B(f₁, f₂) Eq. (142), can be calculated in terms of χ² distributions. For more details see [59].

It is not difficult to examine non-stationarity of the data. One can divide the data into short segments and for each segment calculate the mean, standard deviation and estimate the spectrum. One can then investigate the variation of these quantities from one segment of the data to the other. This simple analysis can be useful in identifying and eliminating bad data. Another quantity to examine is the autocorrelation function of the data. For a stationary process the autocorrelation function should decay to zero. A test to detect certain non-stationarities used for analysis of econometric time series is the Dickey-Fuller test [34]. It models the data by an autoregressive process and it tests whether values of the parameters of the process deviate from those allowed by a stationary model. A robust test for detecting non-stationarity in data from gravitational-wave detectors has been developed by Mohanty [93]. The test involves applying Student’s t-test to Fourier coefficients of segments of the data. Still another block-normal approach has been studied by McNabb et al. [88]. It identifies places in the data stream where the characteristic statistics of the data change. These change points divide the data into blocks in characteristics are stationary.