Coherent Bayesian inference on compact binary inspirals using a network of interferometric gravitational wave detectors

Presented in this paper is a Markov chain Monte Carlo (MCMC) routine for conducting coherent parameter estimation for interferometric gravitational wave observations of an inspiral of binary compact objects using data from multiple detectors. The MCMC technique uses data from several interferometers and infers all nine of the parameters (ignoring spin) associated with the binary system, including the distance to the source, the masses, and the location on the sky. The Metropolis-algorithm utilises advanced MCMC techniques, such as importance resampling and parallel tempering. The data is compared with time-domain inspiral templates that are 2.5 post-Newtonian (PN) in phase and 2.0 PN in amplitude. Our routine could be implemented as part of an inspiral detection pipeline for a world wide network of detectors. Examples are given for simulated signals and data as seen by the LIGO and Virgo detectors operating at their design sensitivity.


I. INTRODUCTION
The era of gravitational wave astronomy is now close upon us as numerous interferometric detectors are operating. The Laser Interferometer Gravitational Wave Observatory (LIGO) [1,2,3] has now reached its target sensitivity, and there is the hope that a detection could come at any time [4]. Around the globe a world-wide network of detectors is coming on-line; Virgo in Italy [5,6,7], GEO in Germany [8,9], and TAMA in Japan [10,11] are operating alongside LIGO in the quest for gravity wave detection. These ground based laser interferometers are sensitive to gravitational radiation in the frequency band from 40 Hz up to 8 kHz.
Coalescing binaries containing neutron stars or black holes promise to be one of the cleanest and most probable sources of detectable radiation [12]. Observation of inspiral events could provide important information on the structure of neutron stars [13,14]. Even cosmological information can be extracted from the observation of inspiral events [15,16,17,18,19]. The characteristics of radiation in the post-Newtonian regime will provide insight into highly non-linear general relativistic effects, such as the observation of the formation of a Kerr black hole as the binary system decays [18,20,21]. The LIGO Scientific Collaboration (LSC) has been actively searching for binary inspiral events [22,23], as well as conducting searches in coincidence with TAMA [24]. Using the data from LIGO's S2 run, it was possible to set an upper limit on the neutron star coalescence rate of less than 50 per year per Milky Way equivalent galaxy [23]. The LSC has also conducted searches for binary inspiral signals from primordial black holes (0.2-1.0 M ⊙ ) in the halo of our galaxy [25], plus more massive black hole systems where component masses are in the 3-20 M ⊙ range [26].
Bayesian inferential methods provide a means to use data from interferometric gravitational wave detectors in order to extract the parameters of a binary inspiral event. Markov chain Monte Carlo (MCMC) methods are a powerful computation technique for parameter estimation within this framework; they are especially useful in applications where the number of parameters is large [27]. Nice descriptions of the positive aspects of a Bayesian analysis of scientific and astrophysical data are provided in [28,29,30]. In previous work we have developed MCMC routines for extracting five parameters associated with a binary inspiral event from data generated by a single interferometric detector [31,32,33]. Our MCMC code took data from a single interferometer, Fourier transformed it into the frequency domain, and then compared the result with 2.0 post-Newtonian (PN) stationary phase templates [34]. One of the new methods that we implement in this current study, presented in this paper, is an MCMC routine that takes time domain interferometric data, and compares it to time domain templates that are 2.5 PN in phase, and 2.0 PN in amplitude; a trivial modification of the code (though not implemented in the study presented here) is to extend the accuracy of the signal waveforms to 3.5 PN in phase and 2.5 PN in amplitude [35,36,37,38,39]. A critical task for a world-wide gravity wave detection network will be to not only detect a binary inspiral signal, but to say where it came from. For this purpose the LSC has developed a coherent binary inspiral search pipeline [40,41,42]. Coherent binary search pipelines and methods are also being developed within the Virgo collaboration [43] and TAMA [44]. Along similar lines, we have developed a coherent MCMC parameter estimation routine, and in this Typeset by REVT E X present paper we describe it and provide results from a test on simulated data. The simulations involve binary neutron star inspirals observed by three well-separated interferometers: the 4 km LIGO detectors at Hanford, WA and Livingston, LA, plus the 3 km Virgo detector in Cascina, Pisa, Italy. The synthetic data for the LIGO and Virgo detectors has Gaussian noise with power spectral densities (PSD) that match their target sensitivities [45,46]. The MCMC code takes data from several interferometers, and estimates the two individual masses, time and phase at coalescence, distance to source, gravity wave polarization angle, angle of inclination of the binary system's orbital plane, and sky position in right ascension and declination. The additional parameters of polarization, inclination angle, and sky location can only be estimated accurately when data from more than one interferometer are considered; they also greatly inflate the parameter space and therefore complicate the analysis. MCMC methods have also been tested in a similar context to recover the nine parameters of a binary black hole coalescence in LISA data [47]. However, the problem setting is different (due to the different instrument, longer observation period, lower frequencies being investigated, and the 2.0 PN phase model), and MCMC techniques were employed rather for optimization than for integration; the MCMC was also only applied on a subset of the parameters, while others were solved for analytically.
The organisation of this paper is as follows. After a brief introduction to the analysis problem we describe our approach alongside more detail on the applied model in Sec. II. Sec. III provides practical directions how we implemented MCMC methods in order to analyze data within the described framework. Sec. IV eventually illustrates results of applications of our code to simulated data. We conclude the paper with a discussion and outlook in Sec. V.

A. Measuring gravitational waves
In an inspiralling binary system, the two companions orbit around their centre of mass with decreasing orbital distance and period, until the system eventually collapses. The gravitational radiation emitted by the system exhibits a 'chirp' form, that is, an oscillation of increasing frequency and amplitude. A laser interferometer is sensitive to space distortions along the directions of its two orthogonal arms, as it monitors the phase difference between the two laser beams that travel along the arms. A gravitational wave is a quadrupole wave that is characterized by its direction of travel, polarization angle, and its two polarization amplitudes. Its effect on a laser interferometer's measurements then is a linear combination of the effects associated with the two polarizations, depending on the orientation of the interferometer with respect to the passing wave. Actual measurements, of course, are also exposed to noise.
Measurements of a binary inspiral's chirp signal by a single interferometer will not be sufficient to infer all of the parameters that determine the signal's waveform and the interferometer's response. Measurements from several separate interferometers will, in general, be required to derive (for example) the wave's direction of travel by matching possible mutual orientations as well as different arrival times of the signals at the different sites. Combining measurements from several interferometers will also enhance sensitivity and signal-to-noise ratio [40].

B. The Bayesian approach
We apply a Bayesian approach to this inferential problem, that is, the term 'probability' is used in a broader sense than in the more common 'frequentist ' interpretation [30,48,49]. Probability calculus here is applied to process and infer states of incomplete information that are reflected by probability distributions, and that are conditional on prior knowledge and/or the data at hand. This allows one to treat unknown parameters as random variables that follow a prior distribution representing the researcher's pre-experimental knowledge and uncertainty. The gain in information through observation of data then follows in a straight-forward fashion through the application of Bayes' theorem, yielding the posterior distribution of the parameters. The posterior distribution, which is essentially the product of the parameters' prior distribution and the likelihood of the data, then poses the basis for inference [50].
Inference through the posterior distribution usually involves the solution of integrals, since one is typically interested in figures such as the marginal (posterior) expectations of individual parameters, marginal (posterior) densities, or (posterior) probabilities of certain events, which are derived from the posterior distribution by integrating over the parameter space. In many cases when analytic integration is not possible, numerical methods are employed, usually (pseudo-) stochastic techniques like Markov chain Monte Carlo (MCMC) methods that simulate random draws from the posterior ditribution, then allowing one to approximate the desired integrals by sample statistics [27,50].

C. Parameters
The waveform that is measured at a certain interferometer depends on the characteristics of the inspiral event as well as the orientation of the source relative to the interferometer. The nine 'global' parameters determining the response of Earth-bound interferometers are: -individual masses (m 1 , m 2 ∈ R + ; m 1 ≤ m 2 ), -luminosity distance (d L ∈ R + ), -inclination angle (ι ∈ [0, π]), -coalescence phase (φ 0 ∈ [0, 2π]), -coalescence time at geocenter (t c ∈ R), -declination (δ ∈ [− π 2 , π 2 ]), -right ascension (α ∈ [0, 2π]) and -polarization (ψ ∈ [0, π]), the latter four of which affect the measurement at the I-th detector in terms of the 'local' parameters -local coalescence time (t These 'local' parameters are derived from the locations/orientations of the source and the individual interferometers with respect to each other. For more specific definitions and conventions see e.g. [51]. In the following we will refer to the two parameter sets as the global parameter vector and the local parameter vector with respect to a specific interferometer I. Not all of the above parameters will usually be of primary interest; especially coalescence phase φ 0 , polarization ψ or inclination ι might be regarded as nuisance parameters.

D. Network likelihood
An interferometer's data output z is assumed to be the sum of the actual signal s(θ), depending on the true parameters θ, and (interferometer-specific) coloured noise. The (real-valued) data z and signal waveform s(θ) enter the likelihood function in terms of their (complexvalued) Fourier transformsz ands(θ), the noise is specified through its power spectral density (PSD) S n . The likelihood function for a specific interferometer I then is (up to a normalising constant) proportional to the following expression [52]. For actual data, discretized and measured over a finite interval of length δ t , it is computed as the sum of squared and normalized differences between the Fourier transforms of the observed signal (z) and the signal template (s(θ)) over the discrete set of Fourier frequencies where i L × ∆ f and i U × ∆ f are the lower and upper bounds of the examined frequency range. Note that, although not labeled as such here, data z, noise spectrum S n , etc. are specific for the different interferometers I. Assuming that noise is independent across different interferometers, the network likelihood then is the product of the individual interferometer likelihoods: (4)

E. Prior specification
The prior information is specified in a straightforward fashion for the 'geometrical' parameters that determine the location and orientation of the inspiral event. A priori, the event is assumed to be equally likely across all possible directions; this leads to a uniform prior for the right ascension α, and a prior density that is proportional to the circumference of the corresponding 'circle of latitude' for the declination δ. Analogously, the prior density of the inclination ι is and the remaining parameters, polarization ψ and coalescence phase φ 0 , have uniform priors. The prior specifications for these parameters may also be regarded as Maximum Entropy choices [49,53].
The coalescence time t c is assumed to be known in advance up to a certain accuracy through preprocessing of the data [22,54,55]; for now we set its prior to be uniform across ±5ms around the true value, which of course is known for simulated data.
The joint prior of the remaining parameters, masses m 1 , m 2 and luminosity distance d L , is set in order to reflect the distribution of parameters given the event has been detected in the first place. Initially, the prior for the two inspiral companions' individual masses is assumed to be uniform between 0.6 and 3.0 M ⊙ (solar masses: M ⊙ ≈ 2×10 30 kg), which effectively covers the range of values expected for neutron stars. The prior for d L is derived from the assumption that inspirals happen uniformly across space, so that P(d L ≤ x) ∝ x 3 . So far, this leads to an improper distribution (that has an infinite integral). But inspiral signals obviously cannot originate from arbitrarily great distances, since at some point their signals become too faint to be detected. We incorporated this restriction by taking into consideration the detection probability D, which we assume to depend on the signal's amplitude, which is roughly proportional to neglecting for simplicity the effects of orientation parameters (A is actually the logarithmic amplitude) [33]. We could have set a threshold amplitude below which inspirals would be assumed undetectable, but favoured a 'smoother' transition that does not explicitly apply zero probability to parts of the parameter space. Instead we model the dependence between signal amplitude and detection probability using a (sigmoidal) logistic function of the form whose parameters a and b are set so that D a,b (x U ) = 1−p and D a,b (x L ) = p, for some upper and lower reference points x U and x L , and some 0 < p < 0.5 (e.g. p := 0.1). So x U denotes the amplitude at which the detection probability reaches 1 − p, and x L is the amplitude at which the probability falls below p. In order to fit D through these points, its parameters are set to So the density of the resulting (proper) prior distribution of individual masses and distance is 0M ⊙ , 50Mpc) and p := 0.1, so an optimally oriented 2.0-2.0 M ⊙ inspiral is assumed to be detectable out to 45 and 50 Mpc with probabilities of 90% and 10%, respectively.
As a 'side effect' of this prior definition, larger masses have a greater prior probability, since inspirals involving large masses may originate from farther distances while low-mass inspirals need to be close in order to be observable at all. This feature is also known as the Malmquist effect ; incorporating it into the prior will compensate for selection bias that would otherwise affect the results [56,57]. The definition of priors, especially for coalescence time t c , individual masses m 1 and m 2 and their relation to the luminosity distance d L , and possibly also for the sky location (δ, α), may be refined at a later stage when e.g. some substantiated knowledge is available about the performance of the upstream detection pipeline, which might provide rough estimates of some of the parameters together with the detection [22,54,55]. For now we aim for simple and general formulations.

A. General
In order to analyze data in terms of the above framework, we implemented an MCMC sampler in C that is supposed to both find the global mode(s) of the posterior distribution and then 'explore' the distribution, i.e. simulate random draws from the posterior. Data is imported from the Frame format using the Frame Library [58]. Prior to the analysis, the data is filtered and downsampled by a factor of 4 [59,60]. (Pseudo-) random number generation within the MCMC sampler was implemented using Randlib [61].
The MCMC sampler writes only every 25th of the drawn samples to a text file, in order to reduce the effects of subsequently correlated samples, and also to keep the data volume at a reasonable level. The MCMC output then is imported into R, a statistical software, for eventual analysis [62].
The marginal densities that are shown in this paper are kernel density estimates [63]. Two-dimensional densities are illustrated by greyscale plots (with darker areas corresponding to greater densities), and in addition the contour line enclosing the most probable region (accumulating 95% probability) is shown.

B. Likelihood implementation
In order to compute the coherent network-likelihood, first the individual interferometer-likelihoods need to be determined. The primary 'ingredients' for the interferometer-likelihood are -the Fourier-transform of the dataz, -the noise spectrum S n , -the 'local' parameter set θ (I) and -the (Fourier-transformed) signal templates(θ (I) ), the first two of which only need to be determined once at the very beginning of the analysis, while the latter two (in general) need to be re-computed for each likelihood evaluation.
For all (discrete) Fourier-transformations we use the freely available FFTW library [64]. The noise spectrum is estimated from a section of data that is disjoint from the actually analyzed data set [65]. In order to minimize undesirable leakage effects, the data is 'windowed' before Fourier-transformation; using a Hann window for spectrum estimation, and a Tukey window for data and template transformations [66].
Internally, along with the noise spectra, data Fourier transforms etc. corresponding to each of the interferometers I, a set of vectors defining the interferometer's location and orientation is stored. This allows to derive the interferometer-specific parameters (local coalescence time t Template waveformss(θ) are generated in the timedomain and then (numerically) Fourier transformed to the frequency domain. Here we used waveform approximations that are 2.5 PN in phase and 2.0 PN in amplitude. The rather complex expressions for these are omitted here and can be found in [51]. We preferred working in the time-domain, since frequency-domain templates might introduce discrepancies because they are exact analytic Fourier transforms of the 'parametric' waveforms, while the actual data is of finite extent and affected by leakage introduced through the numerical discrete Fourier transformation. When using time-domain templates and Fourier transfoming these, the resulting frequency-domain templates are the more accurate match to the Fourier-transformed data. Another advantage of generating templates in the time-domain is the availability of a wider range of signal waveforms; the extension to higher PN appoximations (e.g. 3.5 PN phase and 2.5 PN amplitude [35,36,37,38,39]) or the consideration of additional parameters (e.g. spin effects [70]) would be straightforward to implement.

D. MCMC implementation
In order to enhance the MCMC sampler's performance we applied reparametrisations to some of the parameters. The individual masses (m 1 , m 2 ) are highly correlated in their posterior distribution [18], making sampling from the original parameters extremely difficult. Re-expressing the masses in terms of chirp mass m c and the (symmetric) mass ratio η, where m c = (m 1 m 2 ) 3/5 (m 1 + m 2 ) 1/5 and η = m 1 m 2 (m 1 + m 2 ) 2 , (11) yields a posterior that is much easier to sample from. We then reparameterized the luminosity distance from d L to ln(d L ), which implicitly yields an unbounded parameter space and proposal step widths that are proportional to the current distance d L . Declination δ and inclination ι were transformed to sin(δ) and cos(ι), which leads to uniform prior distributions over the new parameters.
The MCMC algorithm was implemented as a Metropolis-sampler [27,50] that was extended to a parallel tempering algorithm. The idea of tempering is to consider a smoothed ('tempered') version of the actual objective function (here the posterior distribution), or, analogously, do its exploration (optimization, MCMC sampling,...) following 'relaxed' rules. In optimization contexts, the tempering is often faded out over time, the result being a 'simulated annealing' algorithm, which starts off at a high temperature in order to find the global mode amongst other minor modes, but eventually ends up optimizing the original objective function.
Parallel tempering is a special case of the 'Metropoliscoupled MCMC' (MCMCMC) algorithm, in which several 'tempered' chains are run in parallel, each having a different (constant) temperature. So the tempering does not vary over time, but instead is realised across parallel chains, with additional proposals allowing for swapping between chains. Instead of sampling from the regular posterior distribution with density function f (which is essentially the product of prior π and likelihood L: f (θ) ∝ π(θ)L(θ)), the tempered chains sample from a modified distribution where T ≥ 1 denotes the 'temperature', and for which in the extreme cases T = 1 and T → ∞ the tempered distribution f T equals posterior and prior respectively. Chains running at higher temperatures can be considered as sampling from a 'relaxed' or 'widened' posterior which is then used as proposal distribution (through the swapping between chains) for 'cooler' chains, thus improving both convergence and mixing. The draws from the 'coolest' chain with T = 1 are the only ones that are eventually used for posterior inference.
The starting values for the sampler are determined using importance resampling [50]. In a simplified setting of the problem (considering only 5 parameters and one interferometer [33]), this was sufficient to yield reasonable posterior samples that were close enough to the main mode so the sampler then converged reliably and fast. Due to the much larger parameter space and the computationally more expensive likelihood, ensuring convergence through good starting values is not feasible any more. Instead, convergence is now supported through the use of parallel tempering, which enables the sampler to cross gaps between (local) modes and eventually find the posterior's global mode(s).
As proposal distribution for the MCMC sampler we used a multivariate Student's t-distribution with 3 degrees of freedom. In addition to these 'regular' proposals, sometimes draws from the prior are proposed for some parameters in order to enhance convergence, or steps to 'related' parts of the parameter space, like a step from phase φ 0 to φ 0 ± π, which corresponds to an (almost) equally likely parameter combination if the two masses are (almost) equal. Proposals like these are valid as long as a certain symmetry is maintained, i.e. every proposed step is as likely as the reverse step; otherwise one would need to switch to a Metropolis-Hastings sampler that is able to account for asymmetric steps [27,50].

E. Signal-to-noise ratios
The signal-to-noise ratios (SNR) stated in subsequent sections are defined as follows. The interferometerspecific SNR of a certain signal s(θ ⊕ ) received at interferometer I and embedded in noise with spectral density S n is defined as: S n (f ) df (13) [32]. We computed it-in analogy to Eq. (3)-over the same frequency range that is relevant for the likelihood. The network SNR then is defined as: [18].

A. Overview
In the following two sections we will illustrate results generated by our MCMC implementation on simulated data. Firstly, in section IV.B the evidence on parameters in the measured data is illustrated for one simulated signal by looking at the results of an MCMC run in detail. In the following section IV.C the effect of varying signal properties on the evidence in the data is demonstrated and some peculiarities are pointed out. The posterior distributions of parameters are compared across simulated signals that are observed at different distances (and with that, SNRs), but which are otherwise identical. While the individual SNRs (at each interferometer) for each of these examples are similar, a more extreme example with an almost zero SNR at one of the interferometers is considered as well.

B. Recovering the inspiral's parameters
We simulated an inspiral event that is measured at three interferometer sites, namely Hanford (LHO), Livingston (LLO) and Pisa (V). Due to their different noise characteristics, the frequency ranges of the likelihoods were set to 30-1600 Hz for the Virgo observatory, and 40-1600 Hz for the other two LIGO interferometers (cp. Eq. (3)). The amount of data to consider was set to be the 40 seconds before coalescence for Virgo, and 20 s for the other two. This matches the time an inspiral of this kind spends emitting radiation within the above frequency ranges, and would in a realistic search need to be set either with respect to worst-case considerations, or based on prior information supplied by the detection pipeline. The original sampling rates of the data were 20 000 Hz (V) and 16 384 Hz (LHO, LLO), and the signals were superimposed with Gaussian noise matching the corresponding interferometer's design sensitivities [45,46] Six parallel MCMC chains were run within the parallel tempering scheme. With this amount of data the MCMC code generated some 80 samples per minute on a 3.2 GHz Pentium desktop PC, so considering the parallel chains (six) and the thinned output (every 25th), an actual posterior sample is generated every 113 seconds. The first chain of the parallel tempering MCMC sampler converged after some 75 000 iterations, after 'thinning out' of the samples and discarding the burn-in phase, the resulting posterior sample size was 12 500 samples.  Table I lists some numerical posterior estimates. Although correlations between parameters are already greatly reduced through the reparametrisation, some correlation still remains. FIG. 2 illustrates correlations between two such pairs of parameters, one can see that the 'new' mass parameters m c and η are still dependent (though orders of magnitude less than m 1 and m 2 were), and also that the uncertainty in the luminosity distance d L is tied to the uncertainty in the inclination angle ι, since these two parameters can compensate or mimic each other's effect to a certain degree.
Distributions of other variables derived from the parameters could be investigated, their distributions then depending on the joint distribution of the involved parameters [33]. Examples would be the individual masses (m 1 , m 2 ), or the total mass m t = m 1 + m 2 ; the distribution of m t is narrower than those of both m 1 and m 2 , due to their negative correlation (cp. Table I).

C. Results for varying signal characteristics
We performed additional runs with varying settings of the 'true' parameters of the simulated signal. As one would expect, the precision of parameter estimation is proportional to the signal's strength; Table II shows the standard deviations of some of the parameters' posterior distributions. The posterior is narrowest for a close-by inspiral of high masses, and gets wider for both lower mass or greater distance.
These results are in agreement with earlier estimates of the accuracy to be expected from such parameter estimates [18]. The great difference in relative accuracies of parameters related to phase evolution (like chirp mass m c and reduced mass µ = m1m2 m1+m2 = m t η) versus those affecting the signal's amplitude (like distance d L ) is confirmed, and the correlation between m c and µ is verified as well.
At decreasing SNRs, certain parameters cannot be de-termined unambiguously any more. One example is the inclination angle ι, which still has a 'well-behaving' posterior distribution at 10 Mpc distance (see FIG. 2). For a weaker signal originating from 30 Mpc distance, the distribution then turns bimodal (FIG. 3). The 'orientation' of the inclination angle is not clear any more, the result being two roughly equally likely 'mirror image' solutions with P(ι < π 2 ) ≈ 1 2 ≈ P(ι > π 2 ). Note that the two solutions ι and π − ι correspond to opposite orbital directions (clockwise/counterclockwise), as seen from Earth, which might be of minor interest anyway.
The sky location's posterior also exhibits multiple modes for this weaker signal (FIG. 3). This illustrates some potential pitfalls of Maximum-Likelihood (ML) or Maximum-a-Posteriori (MAP) methods; these would advise picking the highest of the several modes, which might  II: Individual and total SNRs for different signals, and some characteristics of the resulting posterior distributions. The accuracy of some of the parameters is illustrated by the posterior standard deviations for (δ, α), tc, dL, mc and µ (percentages refer to the true value). The correlation coefficient for mc and µ shows the (posterior) interdependence between the two parameters. Our results are consistent with those presented in [18].  [71] just be the narrowest one, but not necessarily the most likely. If one then proceeded by extrapolating the curvature at that mode and deriving error bounds from the Fisher Information matrix, the resulting estimates might not only be far off, but also associated with overestimated accuracies.
We also tried MCMC runs with a modified prior setting; we extended the prior for the coalescence time t c from its original range of ±5ms around the true value to ±27ms, allowing for an additional margin of 22ms, which is the time it takes a gravitational wave to travel from Earth's surface to its center. This setting reflects the case where the inspiral detection pipeline received triggers from less than three interferometer sites, so the signal's arrival time at the geocenter could not be estimated to greater accuracy in advance. The MCMC algorithm is still able to find the mode in the enlarged time parameter range, but takes more iterations to converge.
One scenario in which such an approach would be necessary is when the SNR for one of the interferometers is almost zero. In such a case the data from the interferometer under consideration also would not (directly) contribute to the estimation of phase-and frequency-related parameters, but would still carry information about amplitude-related parameters-by implicitly 'ruling out' those parameter combinations that would have resulted in a response at that interferometer. FIG. 4 shows the sky location posteriors for such a signal, a 1.5-2.0 M ⊙ inspiral at 10 Mpc distance, where the SNRs at the three interferometer sites are: LHO: 9.6, LLO: 13.9, V: 0.18 (network: 16.9). Including the data from the third interferometer (with almost zero SNR) into the analysis still yields a much more accurate estimate of the sky location. Table III compares the resulting parameter accuracies of these two settings. The posteriors for sky location (δ, α) and coalescence time t c , which are closely related, are much narrower when the Virgo data is considered in the analysis, while estimates of the rather phase-and frequency-related parameters m c and µ do not gain from the additional information.
On the one hand, not only a high (total) SNR is desirable but also one that is rather 'evenly spread' over different interferometers. On the other hand even a near-zero SNR at one of the interferometers does not make its measurement useless. Inference on different parameters will be affected to different degrees by such an unbalanced SNR arrangement. We have developed a new MCMC program for estimating the nine parameters associated with an inspiral of compact binary objects from the data coming from a network of gravitational wave interferometers. The determination of the sky location of the source is an important consequence of the procedure. Numerous new features have been implemented in this binary inspiral MCMC. The MCMC uses waveform approximations that are 2.5 PN in phase and 2.0 PN in amplitude [51] (and by the time of final submission of this paper a version using waveforms that are 3.5 PN in phase and 2.5 PN in amplitude [35,36,37,38,39] was running as well). The data from multiple interferometers (two or more) are coherently analyzed in order to produce posterior probability distributions for all nine parameters.
Advanced MCMC techniques were implemented in our program in order to maximize the efficiency of converging to the correct parameter values in the large, 9-dimensional, parameter space. The initial parameter values for the sampler were determined using importance resampling [50]. We recently extended (though not with the results presented in this paper) the parallel tempering algorithm to an Evolutionary MCMC algorithm [72]. This MCMC variety implements proposals that are motivated by genetic algorithms [73], and so recombinations of parameter samples from different MCMC chains are used as proposals in order to improve convergence and mixing.
Another current related research effort is the application of a version of this MCMC code to burst waveforms. This problem is by orders of magnitude computationally less expensive, due to the much shorter duration of the signals. But it appears that on the other hand convergence, i.e. finding the main posterior mode in the parameter space, still poses a major problem. The theoretical background of the various potential burst sources is rather vague, so realistic waveforms and sensible specifications of parameterisations and priors also need to be identified. The results of this study on the MCMC parameter estimation of burst signals using the coherent analysis of multi-interferometer data will be presented in a forthcoming publication.
We are also extending the MCMC techniques used in this study to the application of data analysis for LISA detection of binary inspiral signals. While our present program coherently analyzes data from multiple ground based interferometers, we have found it is a straightforward extension of the code so that we can coherently analyze the time delay interferometry data from LISA. These results are also forthcoming.
Presently LIGO is at its target sensitivity. Virgo is fast approaching its design sensitivity. Using the LIGO-Virgo network it will be possible to observe neutron star binary inspirals out to a distance of 35 Mpc [45,46]. A detection of such an inspiral could occur at any time [4]. As displayed in this paper, our MCMC routine is capable of coherently analyzing the data from the multiple interferometers, and then using it to estimate the nine parameters associated with such a signal.