EXPLORING SHORT GAMMA-RAY BURSTS AS GRAVITATIONAL-WAVE STANDARD SIRENS

Two multi-kilometer interferometric gravitational-wave (GW) detectors are presently in operation: LIGO⁴ and Virgo.⁵ They are sensitive to GWs produced by the coalescence of two neutron stars to a distance of roughly 30 Mpc, and to the coalescence of a neutron star with a 10 M_☉ black hole to roughly 60 Mpc. Over the next several years, these detectors will undergo upgrades which are expected to extend their range by a factor ∼10. Most estimates suggest that detectors at advanced sensitivity should measure at least a few, and possibly a few dozen, binary coalescences every year (e.g., Kopparapu et al. 2008; Abadie et al. 2010).

It has long been argued that neutron star–neutron star (NS–NS) and neutron star–black hole (NS–BH) mergers are likely to be accompanied by a gamma-ray burst (Eichler et al. 1989). Recent evidence supports the hypothesis that many short-hard gamma-ray bursts (SHBs) are indeed associated with such mergers (Fox et al. 2005; Nakar et al. 2006; Berger et al. 2007; Perley et al. 2009). This suggests that it may be possible to simultaneously measure a binary coalescence in gamma rays (and associated afterglow emission) and in GWs (Dietz 2009). The combined electromagnetic and gravitational view of these objects will teach us substantially more than what we learn from either data channel alone. Because GWs track a system's global mass and energy dynamics, it has long been known that measuring GWs from a coalescing binary allows us to determine, in the ideal case, "intrinsic" binary properties such as the masses and spins of its members with exquisite accuracy (Finn & Chernoff 1993; Cutler & Flanagan 1994). As we describe in the following subsection, it has also long been appreciated that GWs can determine a system's "extrinsic" properties (Schutz 1986) such as location on the sky and distance to the source. In particular, the amplitude of a binary's GWs directly encodes its luminosity distance. Direct measurement of a coalescing binary could thus be used as a cosmic distance measure: binary inspiral would be a "standard siren" (the GW equivalent of a standard candle, so-called due to the sound-like nature of GWs) whose calibration depends only on the validity of general relativity (Holz & Hughes 2005; Dalal et al. 2006).

Unfortunately, GWs alone do not measure extrinsic parameters as accurately as the intrinsic ones. As we describe in more detail in the following section, GW observation of a binary measures a complicated combination of its distance, its position on the sky, and its orientation, with overall fractional accuracy ∼1/signal-to − noise. As distance is degenerate with these angles, using GWs to measure absolute distance to a source requires a mechanism to break the degeneracy. Associating the GW coalescence waves with a short-hard gamma-ray burst (SHB) is a near-perfect way to break some of these degeneracies.

In this paper, we explore the ability of the near-future advanced LIGO–Virgo detector network to constrain binary parameters (especially distance), when used in conjunction with electromagnetic observations of the same event (such as an associated SHB). We also examine how well these measurements can be improved if planned detectors in Western Australia (AIGO⁶) and in Japan's Kamioka mine (LCGT⁷) are operational. This paper substantially updates and improves upon earlier work (Dalal et al. 2006, hereafter DHHJ06), using a more sophisticated parameter estimation technique. In the next section, we review standard sirens, and in Section 1.3 we briefly summarize DHHJ06. The next subsection describes the organization and background relevant for the rest of the paper.

It has long been recognized that GW inspiral measurements could be used as powerful tools for cosmology. Schutz (1986) first demonstrated this by analyzing how binary coalescences allow a direct measurement of the Hubble constant; Marković (1993) and Finn & Chernoff (1993) subsequently generalized this approach to include other cosmological parameters. More recently, there has been much interest in the measurements enabled when GWs from a merger are accompanied by a counterpart in the electromagnetic spectrum (Bloom et al. 2009; Phinney 2009; Kulkarni & Kasliwal 2009). In this paper, we focus exclusively on GW observations of binaries that have an independent sky position furnished by electromagnetic observations.

We begin by examining GWs from binary inspiral as measured in a single detector. We only present here the lowest order contribution to the waves; in subsequent calculations our results are taken to higher order (see Section 2.1). The leading waveform generated by a source at luminosity distance D_L, corresponding to redshift z, is given by

$$\begin{eqnarray} h_+ &=& \frac{2(1 + z){\cal M}}{D_L}[\pi (1 + z){\cal M}f]^{2/3}(1 + \cos ^2\iota)\cos 2\Phi _N(t), \nonumber \\ h_\times &=& {-}\frac{4(1 + z){\cal M}}{D_L}[\pi (1 + z){\cal M}f]^{2/3}\cos \iota \sin 2\Phi _N(t), \nonumber \\ \Phi _N(t) &=& \Phi _c - \left[\frac{t_c - t}{5(1 + z){\cal M}}\right]^{5/8},\quad f \equiv \frac{1}{\pi }\frac{d\Phi _N}{dt}. \end{eqnarray} \tag{ 1 }$$

Here Φ_N is the lowest order contribution to the orbital phase, f is the GW frequency, and ${\cal M} = m_1^{3/5} m_2^{3/5}/(m_1 + m_2)^{1/5}$ is the binary's "chirp mass," which sets the rate at which f changes. We use units with G = c = 1; handy conversion factors are M_☉ ≡ GM_☉/c² = 1.47 km, and M_☉ ≡ GM_☉/c³ = 4.92 × 10⁻⁶ s. The angle ι describes the inclination of the binary's orbital plane to our line of sight: $\cos \iota = \mathbf {\skew{-4}\hat{L}}\,{\bm \cdot}\, \mathbf {\hat{n}}$, where $\mathbf {\skew{-4}\hat{L}}$ is the unit vector normal to the binary's orbital plane, and $\mathbf {\hat{n}}$ is the unit vector along the line of sight to the binary. The parameters t_c and Φ_c are the time and orbital phase when f diverges in this model. We expect finite size effects to impact the waveform before this divergence is reached.

A given detector measures a linear combination of the polarizations:

$$\begin{equation} h_{\rm meas} = F_+(\theta, \phi, \psi) h_+ + F_\times (\theta, \phi, \psi) h_\times, \end{equation} \tag{ 2 }$$

where θ and ϕ describe the binary's position on the sky, and the "polarization angle" ψ sets the inclination of the components of $\mathbf {\skew{-4}\hat{L}}$ orthogonal to $\mathbf {\hat{n}}$. The angles ι and ψ fully specify the orientation vector $\mathbf {\skew{-4}\hat{L}}$. For a particular detector geometry, the antenna functions F₊ and F_× can be found in Thorne (1987). In Section 2.2, we give a general form for the gravitational waveform without appealing to a specific detector, following the analysis of Cutler & Flanagan 1994 (hereafter abbreviated CF94).

Several features of Equations (1) and (2) are worth commenting upon. First, note that the phase depends on the redshifted chirp mass. Measuring phase thus determines the combination $(1 + z){\cal M}$ (Finn & Chernoff 1993), not ${\cal M}$ or z independently. To understand this, note that ${\cal M}$ controls how fast the frequency evolves: using Equation (1), we find $\dot{f} \propto f^{11/3}{\cal M}^{5/3}$. The chirp mass enters the system's dynamics as a timescale $\tau _c = G{\cal M}/c^3$. For a source at cosmological distance, this timescale is redshifted; the chirp mass we infer is likewise redshifted. Redshift and chirp mass are inextricably degenerate. This remains true even when higher order effects (see, e.g., Blanchet 2006) are taken into account: parameters describing a binary impact its dynamics as timescales which undergo cosmological redshift, so we infer redshifted values for those parameters. GW observations on their own cannot directly determine a source's redshift.

Next, note that the amplitude depends on $(1 + z){\cal M}$, the angles (θ, ϕ, ι, ψ), and the luminosity distance D_L. Measuring the amplitude thus measures a combination of these parameters. By measuring the phase, we measure the redshifted chirp mass sufficiently well that $(1 + z){\cal M}$ essentially decouples from the amplitude. More concretely, matched filtering the data with waveform templates should allow us to determine the phase with fractional accuracy δΦ/Φ ∼ 1/[(signal-to − noise) × (numberofmeasuredcycles)]; $(1 + z){\cal M}$ should be measured with similar fractional accuracy. NS–NS binaries will radiate roughly 10⁴ cycles in the band of advanced LIGO, and NS–BH binaries roughly 10³ cycles, so the accuracy with which phase and redshifted chirp mass can be determined should be exquisite (Finn & Chernoff 1993, CF94).

Although $(1 + z){\cal M}$ decouples from the amplitude, the distance, position, and orientation angles remain highly coupled. To determine source distance we must break the degeneracy that the amplitude's functional form sets on these parameters. One way to break these degeneracies is to measure the waves with multiple detectors. Studies (Sylvestre 2004; Cavalier et al. 2006; Blair et al. 2008; Fairhurst 2009; Wen & Chen 2010) have shown that doing so allows us to determine the source position to within a few degrees in the best cases, giving some information about the source's distance and inclination.

Perhaps the best way to break some of these degeneracies is to measure the event electromagnetically. An EM signature will pin down the event's position far more accurately than GWs alone. The position angles then decouple, much as the redshifted chirp mass decoupled. Using multiple detectors, we can then determine the source's orientation and its distance. This gives us a direct, calibration-free measure of the distance to a cosmic event. The EM signature may also provide us with the event's redshift, directly putting a point on the Hubble diagram. In addition, if modeling or observation gives us evidence for beaming of the SHB emission, this could strongly constrain the source inclination.

Our goal is to assess how well we can determine the luminosity distance D_L to SHBs under the assumption that they are associated with inspiral GWs. We consider both NS–NS and NS–BH mergers as generators of SHBs, and consider several plausible advanced detector networks: the current LIGO/Virgo network, upgraded to advanced sensitivity; LIGO/Virgo plus the proposed Australian AIGO; LIGO/Virgo plus the proposed Japanese LCGT; and LIGO/Virgo plus AIGO plus LCGT.

The engine of our analysis is a probability function that describes how inferred source parameters $\boldsymbol{\theta }$ should be distributed following GW measurement. (Components θ^a of the vector $\boldsymbol{\theta }$ are physical parameters such as a binary's masses, distance, sky position angles, etc.; our particular focus is on D_L.) Consider one detector which measures a datastream s(t), containing noise n(t) and a GW signal $h(t,{{\hat{\boldsymbol\theta }}})$, where ${\hat{\boldsymbol\theta }}$ describes the source's "true" parameters. In the language of Finn (1992), we assume "detection" has already occurred; our goal in this paper is to focus on the complementary problem of "measurement."

As shown by Finn (1992), given a model for our signal $h(t,\boldsymbol{\theta })$, and assuming that the noise statistics are Gaussian, the probability that the parameters $\boldsymbol{\theta }$ describe the data s is

$$\begin{equation} p(\boldsymbol{\theta } | s) = p_0(\boldsymbol{\theta })\exp \left[-\left((h(\boldsymbol{\theta }) - s) |(h(\boldsymbol{\theta }) - s)\right)/2\right]. \end{equation} \tag{ 3 }$$

The inner product (a|b) describes the noise weighted cross-correlation of a(t) with b(t), and is defined precisely below. The distribution $p_0(\boldsymbol{\theta })$ is a prior probability distribution; it encapsulates what we know about our signal prior to measurement. We define ${\tilde{\boldsymbol\theta }}$ to be the parameters that maximize Equation (3).

DHHJ06 did a first pass on the analysis we describe here. They expanded the exponential to second order in the variables $(\boldsymbol{\theta } - {\hat{\boldsymbol\theta }})$; we will henceforth refer to this as the "Gaussian" approximation (cf. Finn 1992):

$$\begin{eqnarray} &&\exp [-(h(\boldsymbol{\theta }) - s | h(\boldsymbol{\theta }) - s)/2]\nonumber\\ &&\quad\simeq \exp \left[-\frac{1}{2}\left(\frac{\partial h}{\partial \theta ^a} \Biggl | \frac{\partial h}{\partial \theta ^b}\right)\delta \theta ^a \delta \theta ^b\right], \end{eqnarray} \tag{ 4 }$$

where $\delta \theta ^a = \theta ^a - \hat{\theta }^a$. In this limit, ${\tilde{\boldsymbol\theta }} = {\hat{\boldsymbol\theta }}$ (at least for uniform priors). The matrix

$$\begin{equation} \Gamma _{ab} \equiv \left(\frac{\partial h}{\partial \theta ^a} \Biggl | \frac{\partial h}{\partial \theta ^b}\right) \end{equation} \tag{ 5 }$$

is the Fisher information matrix. Its inverse Σ^ab is the covariance matrix. Diagonal entries Σ^aa are the variance of parameter θ^a; off-diagonal entries describe correlations.

The Gaussian approximation to Equation (3) is known to be accurate when the signal-to-noise ratio (S/N) is large. However, it is not clear what "large" really means (Vallisneri 2008). Given current binary coalescence rate estimates, it is expected that most events will come from D_L ∼ afew × 100 Mpc. In such cases, we can expect an advanced detector S/N ∼10. It is likely that this value is not high enough for the "large S/N" approximation to be appropriate.

In this analysis we avoid the Gaussian approximation. We instead use Markov Chain Monte Carlo (MCMC) techniques (in particular, the Metropolis–Hastings algorithm) to explore our parameter distributions. A brief description of this technique is given in Section 3, and described in detail in Lewis & Bridle (2002). We find that the Gaussian approximation to Equation (3) is indeed failing in its estimate of extrinsic parameters (though it appears to do well for intrinsic parameters such as mass).

We begin in Section 2 by summarizing how GWs encode the distance to a coalescing binary. We first describe the post-Newtonian (PN) gravitational waveform we use in Section 2.1, and then describe how that wave interacts with a network of detectors in Section 2.2. Our discussion of the network–wave interaction is heavily based on the notation and formalism used in Section 4 of CF94, as well as the analysis of Anderson et al. (2001). Section 2.2 is sufficiently dense that we summarize its major points in Section 2.3 before concluding, in Section 2.4, with a description of the GW detectors which we include in our analysis.

We outline parameter estimation in Section 3. In Section 3.1, we describe in more detail how to construct the probability distributions describing parameter measurement. We then give, in Section 3.2, a brief description of our selection procedure based on S/N detection thresholds. This procedure sets physically motivated priors for some of our parameters. The MCMC technique we use to explore this function is described in Section 3.3. How to appropriately average this distribution to give "noise averaged" results and to compare with previous literature is discussed in Section 3.4.

In Section 4, we discuss the validation of our code. We begin by attempting to reproduce some of the key results on distance measurement presented in CF94. Because of the rather different techniques used by Cutler & Flanagan, we do not expect exact agreement. It is reassuring to find, nonetheless, that we can reconstruct with good accuracy all of the major features of their analysis. We then examine how these results change as we vary the amplitude (moving a fiducial test binary to smaller and larger distances), as we vary the number of detectors in our network, and as we vary the source's inclination.

Our main results are given in Section 5. We consider several different plausible detector networks and examine measurement errors for two "fiducial" binary systems, comprising either two neutron stars (NS–NS) with physical masses of m₁ = m₂ = 1.4 M_☉, or a neutron star and black hole (NS–BH) system with physical masses m₁ = 1.4 M_☉ and m₂ = 10 M_☉. Assuming a constant comoving cosmological density, we distribute potential GW-SHB events on the sky, and select from this distribution using a detection threshold criterion set for the entire GW detector network. We summarize some implications of our results in Section 6. A more in-depth discussion of these implications, particularly with regard to what they imply for cosmological measurements, will be presented in a companion paper.

Throughout this paper, we use units with G = c = 1. We define the shorthand m_z = (1 + z)m for any mass parameter m.

The inspiral and merger of a compact binary's members can be divided into three consecutive phases. The first and longest is a gradual adiabatic inspiral, when the members slowly spiral together due to the radiative loss of orbital energy and angular momentum. PN techniques (an expansion in gravitational potential M/r, or equivalently for bound systems, orbital speed v²) allow a binary's evolution and its emitted GWs to be modeled analytically to high order; see Blanchet (2006) for a review. When the bodies come close together, the PN expansion is no longer valid, and direct numerical calculation is required. Recent breakthroughs in numerical relativity now make it possible to fully model the strong-field, dynamical merger of two bodies into one; see Pretorius (2005), Shibata & Uryū (2006), and Etienne et al. (2008) for discussion. If the end state is a single black hole, the final waves from the system should be described by a ringdown as the black hole settles down to the Kerr solution.

In this work we are concerned solely with the inspiral, and will accordingly use the PN waveform to describe our waves. In particular, we use the so-called restricted PN waveform; following CF94, the inspiral waveform may be written schematically

$$\begin{equation} h(t) = \mathrm{Re}\left(\sum _{x,m} h^x_m(t)e^{im\Phi _{\mathrm{orb}}(t)}\right){.} \end{equation} \tag{ 6 }$$

Here x indicates PN order (h^x is computed to O(v^2x) in orbital speed), m denotes harmonic order (e.g., m = 2 is quadrupole), and Φ_orb(t) = ∫^tΩ(t')dt' is orbital phase (with Ω(t) the orbital angular frequency). The "restricted" waveform neglects all PN amplitude terms beyond the leading one, and considers only the dominant m = 2 contribution to the phase. The phase is computed to high PN order.

Let the unit vector $\hat{\bf n}$ point to a binary on the sky (so that the waves propagate to us along $-\hat{\bf n}$), and let the unit vector $\skew{-4}\hat{\bf L}$ denote the normal along the binary's orbital angular momentum. The waveform is fully described by the two polarizations:

$$\begin{eqnarray} h_+(t) & = & \frac{2{\mathcal M}_z}{D_L} [\pi {\cal M}_z f(t)]^{2/3} [1+(\mathbf {\skew{-4}\hat{L}}\,{\bm \cdot}\, \mathbf {\hat{n}})^2]\cos [\Phi (t)], \nonumber \\ &\equiv & \frac{4{\mathcal M}_z}{D_L} [\pi {\cal M}_z f(t)]^{2/3} {\cal A}_+(\hat{\bf n},\skew{-4}\hat{\bf L})\cos [\Phi (t)]; \end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray} h_{\times }(t) & = & - \frac{4{\mathcal M}_z}{D_L} [\pi {\cal M}_z f(t)]^{2/3} (\mathbf {\skew{-4}\hat{L}}\,{\bm \cdot}\, \mathbf {\hat{n}}) \sin [\Phi (t)]\;, \nonumber \\ &\equiv & \frac{4{\mathcal M}_z}{D_L} [\pi {\cal M}_z f(t)]^{2/3} {\cal A}_\times (\hat{\bf n},\skew{-4}\hat{\bf L})\sin [\Phi (t)]. \end{eqnarray} \tag{ 8 }$$

Equations (7) and (8) are nearly identical to those given in Equation (1); only the phase Φ(t) is different, as described below. ${\cal M}_z$ is the binary's redshifted chirp mass, D_L is its luminosity distance, and we have written the inclination angle cos ι using the vectors $\hat{\bf n}$ and $\skew{-4}\hat{\bf L}$. The functions ${\cal A}_{+,\times }$ compactly gather all dependence on sky position and orientation. In Section 2.2, we discuss how these polarizations interact with our detectors.

In these forms of h₊ and h_×, the phase is computed to 2nd-post-Newtonian (2PN) order (Blanchet et al. 1995):

$$\begin{eqnarray} &&\Phi (t) = 2\pi \int f(t^{\prime })\,dt^{\prime } = 2\pi \int \frac{f}{df/dt}df, \end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray} \frac{df}{dt} &=& \frac{96}{5}\pi ^{8/3} \mathcal {M}_z^{5/3}f^{11/3}\Bigg[1 - \Bigg(\frac{743}{336} + \frac{11}{4}\eta \Bigg)(\pi M_z f)^{2/3} \nonumber \\ & &+ (4\pi)(\pi M_z f) + \Bigg(\frac{34103}{18144} + \frac{13661}{2016}\eta + \frac{59}{18}\eta ^2 \Bigg)\nonumber\\ &&\times(\pi M_z f)^{4/3}\Bigg]. \end{eqnarray} \tag{ 10 }$$

Higher order results for df/dt are now known (Blanchet et al. 2002a, 2002b, 2004), but 2PN order will be adequate for our purposes. Since distance measurements depend on accurate amplitude determination, we do not need a highly refined model of the wave's phase. The rate of sweep is dominantly determined by the chirp mass, but there is an important correction due to η = μ/M = m₁m₂/(m₁ + m₂)², the reduced mass ratio. Note that η is not redshifted; both μ and M (the reduced mass and total mass, respectively) acquire (1 + z) corrections, so their ratio is the same at all z. Accurate measurement of the frequency sweep can thus determine both ${\cal M}_z$ and η (or ${\cal M}_z$ and μ_z).

We will find it useful to work in the frequency domain, using the Fourier transform ${\tilde{h}}(f)$ rather than h(t):

$$\begin{equation} {\tilde{h}}(f) \equiv \int _{-\infty }^{\infty }\, e^{2\pi i f t}h(t)\, dt\;. \end{equation} \tag{ 11 }$$

An approximate result for ${\tilde{h}}(f)$ can be found using stationary phase (Finn & Chernoff 1993), which describes the Fourier transform when f changes slowly:

$$\begin{eqnarray} &&\tilde{h}_+(f) = \sqrt{\frac{5}{96}}\frac{\pi ^{-2/3} {\mathcal M}_z^{5/6}}{D_L}{\cal A}_+ f^{-7/6} e^{i\Psi (f)} \,, \end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray} &&\tilde{h}_\times (f) = \sqrt{\frac{5}{96}}\frac{\pi ^{-2/3} {\mathcal M}_z^{5/6}}{D_L}{\cal A}_\times f^{-7/6} e^{i\Psi (f) - i\pi /2} \,. \end{eqnarray} \tag{ 13 }$$

"Slowly" means that f does not change very much over a single wave period 1/f, so that (df/dt)/f ≪ f. The validity of this approximation for the waveforms we consider, at least until the last moments before merger, has been demonstrated in previous work (Droz et al. 1999). The phase function Ψ(f) in Equations (12) and (13) is given by

$$\begin{eqnarray} \Psi (f) &=& 2\pi f t_c - \Phi _c - \frac{\pi }{4} + \frac{3}{128}(\pi {\mathcal M} f)^{-5/3}\nonumber\\ &&\times \Bigg[1 + \frac{20}{9}\Bigg(\frac{743}{336} + \frac{11}{4}\eta \Bigg)(\pi M_z f)^{2/3} \nonumber \\ & & -16\pi (\pi M_z f) + 10\Bigg(\frac{3058673}{1016064} + \frac{5429}{1008}\eta + \frac{617}{144}\eta ^2 \Bigg)\nonumber\\ &&\times(\pi M_z f)^{4/3}\Bigg]. \end{eqnarray} \tag{ 14 }$$

As in Equation (1), t_c is called the "time of coalescence" and defines the time at which f diverges within the PN framework; Φ_c is similarly the "phase at coalescence." We assume an abrupt and unphysical transition between inspiral and merger at the innermost stable circular orbit (ISCO), $f_{\rm ISCO}=(6 \sqrt{6} \pi M_z)^{-1}$. For NS–NS, f_ISCO occurs at high frequencies where detectors have poor sensitivity. As such, we are confident that this abrupt transition has little impact on our results. For NS–BH, f_ISCO is likely to be in a band with good sensitivity, and better modeling of this transition will be important.

In this analysis we neglect effects which depend on spin. In general relativity, spin drives precessions which can "color" the waveform in important ways, and which can have important observational effects (see, e.g., Vecchio 2004; Lang & Hughes 2006; van der Sluys et al. 2008). These effects are important when the dimensionless spin parameter, a ≡ c|S|/GM², is fairly large. Neutron stars are unlikely to spin fast enough to drive interesting precession during the time that they are in the band of GW detectors. To show this, write the moment of inertia of a neutron star as

$$\begin{equation} I_{\rm NS} = \frac{2}{5}\kappa M_{\rm NS} R_{\rm NS}^2, \end{equation} \tag{ 15 }$$

where M_NS and R_NS are the star's mass and radius, and the parameter κ describes the extent to which its mass is centrally condensed (compared to a uniform sphere). Detailed calculations with different equations of state indicate κ ∼ 0.7–1 (cf. Cook et al. 1994, especially the slowly rotating configurations in their Tables 12, 15, 18, and 21). For a neutron star whose spin period is P_NS, the Kerr parameter is given by

$$\begin{eqnarray} a_{\rm NS} &=& \frac{c}{G}\frac{I_{\rm NS}}{M_{\rm NS}^2} \frac{2\pi }{P_{\rm NS}} \nonumber \\ &\simeq & 0.06 \kappa \left(\frac{R_{\rm NS}}{\rm {12\,km}}\right)^2\left(\frac{1.4\,M_\odot }{M_{\rm NS}}\right) \left(\frac{10\,{\rm ms}}{P_{\rm NS}}\right)\;. \end{eqnarray} \tag{ 16 }$$

As long as the neutron star spin period is longer than ∼10 ms, a_NS is small enough that spin effects can be neglected in our analysis. We should include spin in our models of BH–NS binaries; we leave this to a later analysis. Van der Sluys et al. (2008) included black hole spin effects in an analysis which did not assume a known source position. They found that spin-induced modulations could help GW detectors to localize a source. This and companion works (Raymond et al. 2009; van der Sluys et al. 2009) suggest that, if position is known, spin modulations could improve our ability to measure source inclination and distance.

Our GWs depend on nine parameters: two masses ${\cal M}_z$ and μ_z, two sky position angles (which set $\hat{\bf n}$), two orientation angles (which set $\skew{-4}\hat{\bf L}$), time at coalescence t_c, phase at coalescence Φ_c, and luminosity distance D_L. When the sky position is known, the parameter set is reduced to seven: $\lbrace {\cal M}_z, \mu _z, D_L, t_c, \cos \iota, \psi, \Phi _c \rbrace$.

We now examine how the waves described in Section 2.1 interact with a network of detectors. We begin by introducing a geometric convention, which follows that introduced in CF94 and in Anderson et al. (2001). A source's sky position is given by a unit vector $\hat{\bf n}$ (which points from the center of the Earth to the binary), and its orientation is given by a unit vector $\skew{-4}\hat{\bf L}$ (which points along the binary's orbital angular momentum). We construct a pair of axes which describe the binary's orbital plane:

$$\begin{equation} \hat{\bf X} = \frac{\hat{\bf n}\times \skew{-4}\hat{\bf L}}{|\hat{\bf n}\times \skew{-4}\hat{\bf L}|}\;,\quad \hat{\bf Y} = -\frac{\hat{\bf n}\times \hat{\bf X}}{|\hat{\bf n}\times \hat{\bf X}|}\;. \end{equation} \tag{ 17 }$$

With these axes, we define the polarization basis tensors

$$\begin{eqnarray} &&{\bf e}^+ = \hat{\bf X} \otimes \hat{\bf X} - \hat{\bf Y} \otimes \hat{\bf Y}\;, \end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray} &&{\bf e}^\times = \hat{\bf X} \otimes \hat{\bf Y} + \hat{\bf Y} \otimes \hat{\bf X}\;. \end{eqnarray} \tag{ 19 }$$

The transverse-traceless metric perturbation describing our source's GWs is then

$$\begin{equation} h_{ij} = h_+ e^+_{ij} + h_\times e^\times _{ij}\;. \end{equation} \tag{ 20 }$$

We next characterize the GW detectors. Each detector is an L-shaped interferometer whose arms define two-thirds of an orthonormal triple. Denote by $\hat{\bf x}_a$ and $\hat{\bf y}_a$ the unit vectors along the arms of the ath detector in our network; we call these the x- and y-arms. (The vector $\hat{\bf z}_a = \hat{\bf x}_a \times \hat{\bf y}_a$ points radially from the center of the Earth to the detector's vertex.) These vectors define the response tensor for detector a:

$$\begin{equation} D^{ij}_a = \frac{1}{2}[(\hat{\bf x}_a)^i (\hat{\bf x}_a)^j - (\hat{\bf y}_a)^i (\hat{\bf y}_a)^j]. \end{equation} \tag{ 21 }$$

The response of detector a to a GW is given by

$$\begin{eqnarray} h_a &=& D^{ij}_a h_{ij} \nonumber \\ &\equiv & e^{- 2 \pi i ({\bf n}\,{\bm \cdot}\, {\bf r}_a) f} (F_{a,+}h_+ + F_{a,\times }h_\times), \end{eqnarray} \tag{ 22 }$$

where $\bf {r}_a$ is the position of the detector a and the factor (n · r_a) measures the time of flight between it and the coordinate origin. The second form of Equation (22) shows how the antenna functions introduced in Equation (2) are built from the wave tensor and the response tensor.

Our discussion has so far been frame independent, in that we have defined all vectors and tensors without reference to coordinates. We now introduce a coordinate system for our detectors following Anderson et al. (2001) who in turn use the WGS-84 Earth model (Althouse et al. 1998). The Earth is taken to be an oblate ellipsoid with semi-major axis a = 6.378137 × 10⁶ and semi-minor axis b = 6.356752314 × 10⁶ m. Our coordinates are fixed relative to the center of the Earth. The x-axis (which points along i) pierces the Earth at latitude 0° north, longitude 0° east (normal to the equator at the prime meridian); the y-axis (along j) pierces the Earth at 0° north, 90° east (normal to the equator in the Indian ocean somewhat west of Indonesia); and the z-axis (along k) pierces the Earth at 90° north (the north geographic pole).

A GW source at (θ, ϕ) on the celestial sphere has sky position vector $\hat{\bf n}$:

$$\begin{equation} \hat{\bf n} = \sin \theta \cos \phi {\bf i} + \sin \theta \sin \phi {\bf j} + \cos \theta {\bf k}\;. \end{equation} \tag{ 23 }$$

The polarization angle, ψ, is the angle (measured clockwise about $\hat{\bf n}$) from the orbit's line of nodes to the source's $\hat{\bf X}$-axis. In terms of these angles, the vectors $\hat{\bf X}$ and $\hat{\bf Y}$ are given by (Anderson et al. 2001)

$$\begin{eqnarray} \hat{\bf X} &=& (\sin \phi \cos \psi - \sin \psi \cos \phi \cos \theta) {\bf i}\nonumber\\ &&-(\cos \phi \cos \psi + \sin \psi \sin \phi \cos \theta){\bf j}\nonumber\\ &&+ \sin \psi \sin \theta {\bf k}, \end{eqnarray} \tag{ 24 }$$

$$\begin{eqnarray} \hat{\bf Y} &=& ({-}\sin \phi \sin \psi - \cos \psi \cos \phi \cos \theta) {\bf i}\nonumber\\ &&+ (\cos \phi \sin \psi - \cos \psi \sin \phi \cos \theta){\bf j}\nonumber\\ &&+ \cos \psi \sin \theta {\bf k}. \end{eqnarray} \tag{ 25 }$$

The angle ϕ is related to right ascension α by α = ϕ + GMST (where GMST is the Greenwich mean sidereal time at which the signal arrives), and θ is related to declination δ by δ = π/2 − θ (cf. Anderson et al. 2001, Appendix B). Combining Equations (24) and (25) with Equations (18)–(20) allows us to write h_ij for a source in coordinates adapted to this problem.

We now similarly describe our detectors using convenient coordinates. Detector a is at east longitude λ_a and north latitude φ_a (not to be confused with sky position angle ϕ). The unit vectors pointing east, north, and up for this detector are

$$\begin{eqnarray} &&{\bf e}^{\rm E}_a = {-}\sin \lambda _a{\bf i} + \cos \lambda _a{\bf j}, \end{eqnarray} \tag{ 26 }$$

$$\begin{eqnarray} &&{\bf e}^{\rm N}_a = -\sin \varphi _a\cos \lambda _a{\bf i} - \sin \varphi _a\sin \lambda _a{\bf j} + \cos \varphi _a{\bf k}\;, \end{eqnarray} \tag{ 27 }$$

$$\begin{eqnarray} &&{\bf e}^{\rm U}_a = \cos \varphi _a\cos \lambda _a{\bf i} + \cos \varphi _a\sin \lambda _a{\bf j} - \cos \varphi _a{\bf k}\;. \end{eqnarray} \tag{ 28 }$$

The x-arm of detector a is oriented at angle ϒ_a north of east, while its y-arm is at angle ϒ_a + π/2. Thanks to Earth's oblateness, the x- and y-arms are tilted at angles ω^x,y_a to the vertical. The unit vectors $\hat{\bf x}_a$, $\hat{\bf y}_a$ can thus be written as

$$\begin{eqnarray} &&\hat{\bf x}_a = \cos \omega ^x_a \cos \Upsilon _a{\bf e}^{\rm E}_a + \cos \omega ^x_a\sin \Upsilon _a{\bf e}^{\rm N}_a + \sin \omega ^x_a{\bf e}^{\rm U}, \end{eqnarray} \tag{ 29 }$$

$$\begin{eqnarray} &&\hat{\bf y}_a = {-}\cos \omega ^y_a\sin \Upsilon _a{\bf e}^{\rm E}_a + \cos \omega ^y_a\cos \Upsilon _a{\bf e}^{\rm N}_a + \sin \omega ^y_a{\bf e}^{\rm U}. \end{eqnarray} \tag{ 30 }$$

Combining Equations (29) and (30) with Equation (21) allows us to write the response tensor for each detector in our network.

Section 2.2 is sufficiently dense that a brief summary may clarify its key features, particularly with respect to the quantities we hope to measure. From Equation (22), we find that each detector in our network measures a weighted sum of the two GW polarizations h₊ and h_×. Following Cutler (1998), we can rewrite the waveform detector a measures as

$$\begin{equation} h_a = \frac{4{\cal M}_z}{D_L}{\cal A}_p\left[\pi {\cal M}_z f(t)\right]^{2/3} \cos [\Phi (t) + \Phi _p], \end{equation} \tag{ 31 }$$

where we have introduced detector a's "polarization amplitude"

$$\begin{equation} {\cal A}_p = \sqrt{ (F_{a,+}{\cal A}_+)^2 + (F_{a,\times }{\cal A}_\times)^2}, \end{equation} \tag{ 32 }$$

and its "polarization phase"

$$\begin{equation} \tan \Phi _p = \frac{F_{a,\times }{\cal A}_\times }{F_{a,+}{\cal A}_+}\;. \end{equation} \tag{ 33 }$$

The intrinsic GW phase, Φ(t), is a strong function of the redshifted chirp mass, ${\cal M}_z$, the redshifted reduced mass, μ_z, the time of coalescence, t_c, and the phase at coalescence, Φ_c. Measuring the phase determines these four quantities, typically with very good accuracy.

Consider for a moment measurements by a single detector. The polarization amplitude and phase depend on the binary's sky position, (θ, ϕ) or $\hat{\bf n}$, and orientation, (ψ, ι) or $\skew{-4}\hat{\bf L}$. (They also depend on detector position, (λ_a, φ_a), orientation, ϒ_a, and tilt, (ω^x_a, ω^y_a). These angles are known and fixed, so we ignore them in this discussion.) If the angles (θ, ϕ, ψ, ι) are not known, a single detector cannot separate them, nor can it separate the distance D_L.

Multiple detectors can, at least in principle, separately determine these parameters. Each detector measures its own amplitude and polarization phase. Combining their outputs, we can fit to the unknown angles and the distance. Various works have analyzed how well this can be done assuming that the position and orientation are completely unknown (Sylvestre 2004; Cavalier et al. 2006; Blair et al. 2008). Van der Sluys et al. (2008) performed such an analysis for measurements of NS–BH binaries, including the effect of orbital precession induced by the black hole. This precession effectively makes the angles ι and ψ time dependent, breaking the degeneracy among these angles and D_L.

In what follows, we assume that an electromagnetic identification pins down the angles (θ, ϕ), so that they do not need to be determined from the GW data. We then face the substantially less challenging problem of determining ψ, ι, and D_L. We will also examine the impact of a constraint on the inclination, ι. Long bursts are believed to be strongly collimated, emitting into jets with opening angles of just a few degrees. Less is known about the collimation of SHBs, but it is plausible that their emission may be primarily along a preferred axis (presumably the progenitor binary's orbital angular momentum axis).

Here we briefly summarize the properties of the GW detectors that we consider.

LIGO. The Laser Interferometer Gravitational-wave Observatory consists of two 4 km interferometers located in Hanford, Washington (US) and Livingston, Louisiana (US). These instruments have achieved their initial sensitivity goals. An upgrade to "advanced" configuration is expected to be completed around 2014, with tuning for best sensitivity to be undertaken in the years following.⁸ We show the anticipated noise limits from fundamental noise sources in Figure 1 for a broadband tuning (Harry & LIGO Scientific Collaboration 2010). This spectrum is expected to be dominated by quantum sensing noise above a cutoff at f < 10 Hz, with a contribution from thermal noise in the test mass coatings in the band from 30 to 200 Hz.

Zoom In Zoom Out Reset image size

Virgo. The Virgo detector (Acernese & Virgo Scientific Collaboration 2008) near Pisa, Italy has slightly shorter arms than LIGO (3 km), but should achieve similar advanced sensitivity on roughly the same timescale as the LIGO detectors.⁹ For simplicity, we will take Virgo's sensitivity to be the same as LIGO's.

Our baseline detector network consists of the LIGO Hanford and Livingston sites, and Virgo; these are instruments which are running today, and will be upgraded over the next decade. We also examine the impact of adding two proposed interferometers to this network.

AIGO. The Australian International Gravitational Observatory (Barriga et al. 2010) is a proposed multi-kilometer interferometer that would be located in Gingin, Western Australia. AIGO's proposed site in Western Australia is particularly favorable due to low seismic and human activity.

LCGT. The Large-scale Cryogenic Gravitational-wave Telescope (Kuroda & LCGT Collaboration 2010) is a proposed multi-kilometer interferometer that would be located in the Kamioka observatory, 1 km underground. This location takes advantage of the fact that local ground motions tend to decay rapidly as we move away from Earth's surface. They also plan to use cryogenic cooling to reduce thermal noise.

As with Virgo, we will take the sensitivity of AIGO and LCGT to be the same as LIGO for our analysis. Table 1 gives the location and orientation of these detectors, needed to compute each detector's response function. It is worth mentioning that more advanced detectors are in the early planning stages. Particularly noteworthy is the European proposal for the "Einstein Telescope," currently undergoing design studies. It is being designed to study binary coalescence to high redshift (z ≳ 5) (Sathyaprakash et al. 2009).

Table 1. GW Detectors (Positions and Orientations)

Detector	East Long. λ	North Lat. φ	Orientation ϒ	x-arm Tilt ωx	y-arm Tilt ωy
LIGO-Han	−1194	465	126°	(−6.20 × 10−4)°	(1.25 × 10−5)°
LIGO-Liv	−908	306	198°	(−3.12 × 10−4)°	(−6.11 × 10−4)°
Virgo	105	436	70°	00	00
AIGO	1157	−314	0°	00	00
LCGT	1373	364	25°	00	00

Download table as:  ASCIITypeset image

We now give a brief summary of the parameter estimation formalism we use. Further details can be found in Finn (1992), Królak et al. (1993), and CF94.

Assuming detection has occurred, the datastream of detector a, s_a(t), has two contributions: the true GW signal $h_a(t;{{\hat{\boldsymbol\theta }}})$ (constructed by contracting the GW tensor h_ij with detector a's response tensor D^ij_a; cf. Section 2.2), and a realization of detector noise n_a(t),

$$\begin{equation} s_a(t) = h_a(t; {{\hat{\boldsymbol\theta }}}) + n_a(t)\;. \end{equation} \tag{ 34 }$$

The incident GW strain depends on (unknown) true parameters ${{\hat{\boldsymbol\theta }}}$. As in Section 1.3, ${\hat{\boldsymbol\theta }}$ is a vector whose components are binary parameters. Below we use a vector s whose components s_a are the datastreams of each detector. Likewise, h and n are vectors whose components are the GW and noise content of each detector.

We assume the noise to be stationary, zero mean, and Gaussian. This lets us categorize it using the spectral density as follows. First, define the noise correlation matrix:

$$\begin{eqnarray} C_n(\tau)_{ab} &=& \langle n_a(t + \tau) n_b(t) \rangle - \langle n_a(t + \tau) \rangle \, \langle n_b(t) \rangle \nonumber \\ &=& \langle n_a(t + \tau) n_b(t) \rangle \;, \end{eqnarray} \tag{ 35 }$$

where the angle brackets are ensemble averages over noise realizations, and the zero mean assumption gives us the simplified form on the second line. For a = b, this is the auto-correlation of detector a's noise; otherwise, it describes the correlation between detectors a and b. The (one-sided) power spectral density matrix is the Fourier transform of this:

$$\begin{equation} S_n(f)_{ab} = 2 \int _{-\infty }^{\infty } d \tau \, e^{2 \pi i f \tau } C_n(\tau)_{ab}\;. \end{equation} \tag{ 36 }$$

This is defined for f > 0 only. For a = b, it is the spectral density of noise power in detector a; for a ≠ b, it again describes correlations between detectors. From these definitions, one can show that

$$\begin{equation} \langle {\tilde{n}}_a(f) \, {\tilde{n}}_b(f^\prime)^* \rangle = {1 \over 2} \delta (f - f^\prime) S_n(f)_{ab}. \end{equation} \tag{ 37 }$$

For Gaussian noise, this statistic completely characterizes our detector noise. No real detector is completely Gaussian, but by using multiple, widely separated detectors non-Gaussian events can be rejected. For this analysis, we assume that the detectors' noises are uncorrelated such that Equation (37) becomes

$$\begin{equation} \langle {\tilde{n}}_a(f) \, {\tilde{n}}_b(f^\prime)^* \rangle = {1 \over 2} \delta _{ab} \delta (f - f^\prime) S_n(f)_a. \end{equation} \tag{ 38 }$$

Finally, for simplicity we assume that S_n(f)_a has the universal shape S_n(f) projected for advanced LIGO, as shown in Figure 1.

Many of our assumptions are idealized (Gaussian noise; identical noise spectra; no correlated noise between interferometers), and will certainly not be achieved in practice. These idealizations greatly simplify our analysis, however, and are a useful baseline. It would be useful to revisit these assumptions and understand the quantitative impact that they have on our analysis, but we do not expect a major qualitative change in our conclusions.

The central quantity of interest in parameter estimation is the posterior probability distribution function (PDF) for ${\boldsymbol{\theta }}$ given detector output s, which is defined as

$$\begin{equation} p({\boldsymbol{\theta }} \, | \, {\bf s}) = {\cal N} \, p^{(0)} ({\boldsymbol{\theta }}) {\cal L}_{\rm TOT} ({\bf s} \, | \,{\boldsymbol{\theta }}) \,. \end{equation} \tag{ 39 }$$

${\cal N}$ is a normalization constant, $p^{(0)}({\boldsymbol{\theta }})$ is the PDF that represents the prior probability that a measured GW is described by the parameters $\boldsymbol{\theta }$, and ${\cal L}_{\rm TOT} (\bf {s} \, | \, {\boldsymbol{\theta }})$ is the total likelihood function (e.g., MacKay 2003). The likelihood function measures the relative conditional probability of observing a particular data set $\bf {s}$ given a measured signal h depending on some unknown set of parameters $\boldsymbol{\theta }$ and given noise n. Because we assume that the noise is independent and uncorrelated at each detector site, we may take the total likelihood function to be the product of the individual likelihoods at each detector:

$$\begin{equation} {\cal L}_{\rm TOT} ({\bf s} \, | \,{\boldsymbol{\theta }}) = \Pi _{a} {\cal L}_a (s_a \, | \,{\boldsymbol{\theta }})\;, \end{equation} \tag{ 40 }$$

where ${\cal L}_a$, the likelihood for detector a, is given by (Finn 1992)

$$\begin{equation} {\cal L}_a \, (s \, | \,{\boldsymbol{\theta }}) = \, e^{ - (h_a({\boldsymbol{\theta }}) - s_a \, \big | \, h_a({\boldsymbol{\theta }}) - s_a)/2 } \,. \end{equation} \tag{ 41 }$$

The inner product (...|...) on the vector space of signals is defined as

$$\begin{equation} (g|h) = 2 \int _0^{\infty } df \frac{\tilde{g}^*(f)\tilde{h}(f) + \tilde{g}(f)\tilde{h}^*(f)}{S_n(f)} \,. \end{equation} \tag{ 42 }$$

This definition means that the probability of the noise n(t) taking some realization n₀(t) is

$$\begin{equation} p(n = n_0) \, \propto \, e^{- \left(n_0 | n_0 \right) / 2 }. \end{equation} \tag{ 43 }$$

For clarity, we distinguish between various definitions of S/N. The true S/N at detector a, associated with a given instance of noise for a measurement at a particular detector, is defined as (CF94)

$$\begin{eqnarray} \left({S\over N}\right)_{a, {\rm true}} & = & { \left(h_a \, | \, s_a \right) \over \sqrt{ \left(h_a \, | \, h_a\right) } }\;. \end{eqnarray} \tag{ 44 }$$

This is a random variable with Gaussian PDF of unit variance. For an ensemble of realizations of the detector noise n_a, the average S/N at detector a is given by

$$\begin{equation} \left({S\over N}\right)_{a, {\rm ave}} = {{(h_a | h_a)}\over { {\rm rms}\ (h_a|n_a)}} = (h_a|h_a)^{1/2}. \end{equation} \tag{ 45 }$$

Consequently, we can define the combined true and average S/Ns of a coherent network of detectors:

$$\begin{eqnarray} \left({S\over N}\right)_{{\rm true}} & = & \sqrt{\sum _a \left({S\over N}\right)^2_{a, {\rm true}}}\ \, \end{eqnarray} \tag{ 46 }$$

and

$$\begin{eqnarray} \left({S\over N}\right)_{{\rm ave}} & = & \sqrt{\sum _a \left({S\over N}\right)^2_{a, {\rm ave}}}. \end{eqnarray} \tag{ 47 }$$

Estimating the parameter set ${\boldsymbol{\theta }}$ is often done using a "maximum-likelihood" method following either a Bayesian (Loredo 1989; Finn 1992; CF94, Poisson & Will 1995) or frequentist point of view (Królak et al. 1993, CF94). We do not attempt to review these philosophies, and instead refer to Appendix A2 of CF94 for detailed discussion. It is worth noting that, in the GW literature, the "maximum likelihood" or "maximum a posterior" are often interchangeably referred to as "best-fit" parameters. The maximum a posterior is the parameter set ${\tilde{\boldsymbol\theta }}_{\rm MAP}$ which maximizes the full posterior probability, Equation (39); likewise, the maximum likelihood is the parameter set ${\tilde{\boldsymbol\theta }}_{\rm ML}$ which maximizes the likelihood function, Equation (40).

Following the approach advocated by CF94, we introduce the Bayes estimator ${\tilde{\theta }}_{\rm BAYES}^i({\bf s})$,

$$\begin{equation} {\tilde{\theta }}_{\rm BAYES}^i({\bf s}) \equiv \int {\theta }^i\, p(\boldsymbol{\theta } \, | \, {\bf s}) d\boldsymbol{\theta }\;. \end{equation} \tag{ 48 }$$

The integral is performed over the whole parameter set $\boldsymbol{\theta }$; $d\boldsymbol{\theta } = d\theta ^1d\theta ^2\dots d\theta ^n$. Similarly, we define the rms measurement errors Σ^ij_BAYES

$$\begin{equation} \Sigma _{\rm BAYES}^{ij} = \int \big({\theta }^i - {\tilde{\theta }}^i_{\rm BAYES}\big) \, ({\theta }^j - {\tilde{\theta }}^j_{\rm BAYES}\big) \, p(\boldsymbol{ \theta } \, | \, {\bf s}) d\boldsymbol{\theta }. \end{equation} \tag{ 49 }$$

To understand the meaning of ${\tilde{\theta }}_{\rm BAYES}^i({\bf s})$, consider a single detector which records an arbitrarily large ensemble of signals. This ensemble will contain a sub-ensemble in which the various s(t) are identical to one another. Each member of the sub-ensemble corresponds to GW signals with different true parameters ${\hat{\boldsymbol\theta }}$, but have noise realizations n(t) that conspire to produce the same s(t). In this case, ${\tilde{\theta }}_{\rm BAYES}^i({\bf s})$ is the expectation of θⁱ averaged over the sub-ensemble. The principal disadvantage of the Bayes estimator is the computational cost to evaluate the multi-dimensional integrals in Equations (48) and (49).

For large S/N it can be shown that the estimators ${\tilde{\boldsymbol\theta }}_{\rm ML}$, ${\tilde{\boldsymbol\theta }}_{\rm MAP}$, and ${\tilde{\boldsymbol\theta }}_{\rm BAYES}$ agree with one another (CF94), and that Equation (39) is well described by a Gaussian form (cf. Equation (4)). However, as illustrated in Section IVD of CF94, effects due to prior information and which scale nonlinearly with 1/S/N contribute significantly at low S/N. The Gaussian approximation then tends to underestimate measurement errors by missing tails or multimodal structure in posterior distributions.

We emphasize that in this analysis we do not consider systematic errors that occur due to limitations in our source model or to gravitational lensing effects. A framework for analyzing systematic errors in GW measurements has recently been presented by Cutler & Vallisneri (2007). An important follow-on to this work will be to estimate systematic effects and determine whether they significantly change our conclusions.

We now describe how we generate a sample of detectable GW-SHB events. We assume a constant comoving density (Peebles 1993; Hogg 1999) of GW-SHB events, in a ΛCDM universe with H₀ = 70.5 kms^-1Mpc^-1, Ω_Λ = 0.726, and Ω_m = 0.2732 (Komatsu et al. 2009). We distribute 10⁶ binaries uniformly in volume with random sky positions and orientations to redshift z = 1 (D_L ≃ 6.6 Gpc). We then compute the average S/N, Equation (45), for each binary at each detector, and use Equation (47) to compute the average total S/N for each network we consider. We assume prior knowledge of the merger time (since we have assumed that the inspiral is correlated with an SHB), so we set a threshold S/N for the total detector network, S/N_total = 7.5 (see discussion in DHHJ06). This is somewhat reduced from the threshold we would set in the absence of a counterpart, since prior knowledge of merger time and source position reduces the number of search templates we need by a factor ∼10⁵ (Kochanek & Piran 1993; Owen 1996). Using the average S/N to set our threshold introduces a slight error into our analysis, since the true S/N will differ from the average. Some events which we identify as above threshold could be moved below threshold due to a measurement's particular noise realization. However, some sub-threshold events will likewise be moved above threshold, and the net effect is not expected to be significant.

Our threshold selects detectable GW-SHB events for each detector network. We define "total detected binaries" to mean binaries which are detected by a network of all five detectors—both LIGO sites, Virgo, AIGO, and LCGT. Including AIGO and LCGT substantially improves the number detected, as compared to just using the two LIGO detectors and Virgo. Assuming that all binary orientations are equally likely given an SHB (i.e., no beaming), we find that a LIGO–Virgo network detects 50% of the total detected binaries; LIGO–Virgo–AIGO detects 74% of the total; and LIGO–Virgo–LCGT detects 72% of the total. Figure 2 shows the sky distribution of detected binaries for various detector combinations. Networks which include LCGT tend to have rather uniform sky coverage. Those with AIGO cover the quadrants cos θ>0, ϕ>π and cos θ < 0, ϕ < π particularly well.

Zoom In Zoom Out Reset image size

Our selection method implicitly sets a prior distribution on our parameters. For example, the thresholding procedure results in a significant bias in detected events toward face-on binaries, with $\mathbf {\skew{-4}\hat{L}}\,{\bm \cdot}\, \mathbf {\hat{n}} \rightarrow \pm 1$. Figure 3 shows the distribution of detectable NS–NS binaries for the parameters (cos ι, D_L). Since we use an unrealistic mass distribution (1.4 M_☉–1.4 M_☉ NS–NS and 1.4 M_☉–10 M_☉ NS–BH binaries), instead of a more astrophysically realistic distribution, the implicit mass prior is uninteresting. Figure 4 shows the average total S/N versus the true D_L of our sample of detectable NS–NS and NS–BH binaries for our "full" network (LIGO, Virgo, AIGO, LCGT). Very few detected binaries have S/N above 30 for NS–NS, and above 70 for NS–BH. It is interesting to note the different detectable ranges between the two populations: NS–BH binaries are detectable to over twice the distance of NS–NS binaries.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We are also interested in seeing the impact that prior knowledge of SHB collimation may have on our ability to measure these events. To date there exist only two tentative observations which suggest that SHBs may be collimated (Grupe et al. 2006; Burrows et al. 2006; Soderberg et al. 2006); we therefore present results for moderate collimation and for isotropic SHB emission. To obtain a sample of beamed SHBs, we assume that the burst emission is collimated along the orbital angular momentum axis, where baryon loading is minimized. Following DHHJ06, we use a distribution for cos ι ≡ v of dP/dv ∝ exp [ − (1 − v)²/2σ²_v], with σ_v = 0.05. This corresponds to a beamed population with 68% of its distribution having an opening jet angle within roughly 25°. We construct a beamed subsample by selecting events from the total sample of detected events such that the final distribution in inclination angle follows dP/dv. Joint measurements of SHBs and GW-driven inspirals should enable us to constrain beaming angles by comparing the measured rates for these two populations.

The principal disadvantage of the Bayes estimators $\tilde{\theta }^i_{\rm BAYES}$ and Σ^ij_BAYES is the high computational cost of evaluating the multi-dimensional integrals which define them, Equations (48) and (49). To get around this problem, we use MCMC methods to explore the PDFs describing the seven parameters $\lbrace {\cal M}_c, \mu, D_L, \cos \iota, \psi, t_c, \Phi _c \rbrace$. MCMC methods are widely used in diverse astrophysical applications, ranging from high-precision cosmology (e.g., Dunkley et al. 2009; Sievers et al. 2009) to extra-solar planet studies (e.g., Ford 2005; Winn et al. 2007). They have seen increased use in GW measurement and parameter estimation studies in recent years (e.g., Stroeer et al. 2006; Wickham et al. 2006; Cornish & Porter 2007; Porter & Cornish 2008; Röver et al. 2007; van der Sluys et al. 2008).

MCMC generates a random sequence of parameter states that sample the posterior distribution, $p(\boldsymbol{\theta } | \mathbf {s})$. Let the nth sample in the sequence be $\boldsymbol{\theta }^{(n)}$. If one draws a total of N random samples, Equations (48) and (49) can then be approximated as sample averages:

$$\begin{eqnarray} &&{\tilde{\theta }}^i_{\rm BAYES} \simeq \frac{1}{N} \sum _{n = 1}^N (\theta ^i)^{(n)}\;, \end{eqnarray} \tag{ 50 }$$

$$\begin{eqnarray} &&\Sigma ^{ij}_{\rm BAYES} \simeq \frac{1}{N} \sum _{n = 1}^N \big(\tilde{\theta }^i_{\rm BAYES} - (\theta ^i)^{(n)}\big) \big(\tilde{\theta }^j_{\rm BAYES} - (\theta ^j)^{(n)}\big).\quad \end{eqnarray} \tag{ 51 }$$

The key to making this technique work is drawing a sequence that represents the posterior PDF. We use the Metropolis–Hastings algorithm to do this (Metropolis et al. 1953; Hastings 1970); see Neal (1993), Gilks et al. (1996), MacKay (2003), and Christensen et al. (2004) for in-depth discussion. The MCMC algorithm we use is based on a generic version of CosmoMC,¹⁰ described in Lewis & Bridle (2002).

Appropriate priors are crucial to any MCMC analysis. We take the prior distributions in chirp mass ${\cal M}_z$, reduced mass μ_z, polarization angle ψ, coalescence time t_c, and coalescence phase Φ_c to be flat over the region of sample space where the binary is detectable according to our selection procedure. More specifically, we choose the following.

• 1.  
  $p^{(0)}({\cal M}_z) = {\rm constant}$ over the range [1 M_☉, 2 M_☉] for NS–NS and over the range [2.5 M_☉, 4.9 M_☉] for NS–BH. (The true chirp masses in the binaries' rest frames are 1.2 M_☉ for NS–NS and 3.0 M_☉ for NS–BH.)
• 2.  
  p⁽⁰⁾(μ_z) = constant over the range [0.3 M_☉, 2 M_☉] for NS–NS and over the range [0.5 M_☉, 3.5 M_☉] for NS–BH. (The true reduced masses in the binaries' rest frames are 0.7 M_☉ for NS–NS and 1.2 M_☉ for NS–BH.)
• 3.  
  p⁽⁰⁾(ψ) = constant over the range [0, π].
• 4.  
  p⁽⁰⁾(t_c) = constant over the range [ − 100 s, 100 s]. Since we assume that t_c is close to the time of the SHB event, it is essentially the time offset between the system's final GWs and its SHB photons. We find that the range in t_c we choose is almost irrelevant, as long as the prior is flat and includes the true value. No matter how broad we choose the prior in t_c, our posterior PDF ends up narrowly peaked around $\hat{t}_c$.
• 5.  
  p⁽⁰⁾(Φ_c) = constant over the range [0, 2π].

The prior distribution for D_L is inferred by taking the density of SHBs to be uniform per unit comoving volume over the luminosity distance range [0, 2 Gpc] for NS–NS binaries, and over the range [0, 5 Gpc] for NS–BH binaries. For our sample with isotropic inclination distribution, we put p⁽⁰⁾(cos ι) = constant over the range [ − 1, 1]. When we assume SHB collimation, our prior in cos ι ≡ v is the same as the one that we used in our selection procedure discussed previously:

$$\begin{equation} \frac{dp^{(0)}}{dv}(v) \propto e^{-(1 -v)^2/2\sigma _v^2}\;, \end{equation} \tag{ 52 }$$

with σ_v = 0.05.

We then map out full distributions for each of our seven parameters, assessing the mean values (Equation (48)) and the standard deviations (Equation (49)). We generate four chains which run in parallel on the CITA "Sunnyvale" Cluster. Each chain runs for a maximum of 10⁷ steps; we find that the mean and median number of steps are ∼10⁵ and ∼10⁴, respectively. Each evaluation of the likelihood function takes ∼0.3 s. We use the first 30% of a chain's sample states for "burn in," and discard that data. Our chains start at random offset parameter values, drawn from Gaussians centered on the true parameter values. We assess convergence by testing whether the multiple chains have produced consistent parameter distributions. Following standard practice, we use the Gelman–Rubin convergence criterion, defining a sequence as "converged" if the statistic R < 1.1 on the last half of our samples; see Gelman & Rubin (1992) for more details. We use convergence as our stopping criterion. Each simulation for every binary runs for 1–48 hr; the mean and median runtime are 8 hr and 3 hr, respectively.

Central to the procedure outlined above is the use of the datastream ${\bf s} = {\bf h}(\boldsymbol{\theta }) + {\bf n}$ which enters the likelihood function ${\cal L}_{\rm TOT}({\bf s} | \boldsymbol{\theta })$. The resulting posterior PDF, and the parameters one infers, thus depends on the noise n which one uses. One may want to evaluate statistics that are in a well-defined sense "typical" given the average noise properties, rather than depending on a particular noise instance. Such averaging is appropriate, for example, when forecasting how well an instrument should be able to measure the properties of a source or a process. We have also found it is necessary to average when trying to compare our MCMC code's output with previous work.

As derived below, the averaged posterior PDF takes a remarkably simple form: it is the "usual" posterior PDF, Equation (39) with the noise n set to zero. This does not mean that one ignores noise when constructing the averaged PDF; one still relates the signal amplitude to typical noise by the average S/N, Equation (45). As such, the averaged statistics will show an improvement in measurement accuracy as S/N is increased.

To develop a useful notion of averaged posterior PDF, consider the hypothetical (and wholly unrealistic) case in which we measure a signal using M different noise realizations for the same event. The joint likelihood for these measurements is

$$\begin{equation} {\cal L}^{\rm joint}_{\rm TOT}({\bf s}_1, {\bf s}_2, \ldots {\bf s}_M | \boldsymbol{\theta }) = \prod _{i=1}^M {\cal L}_{\rm TOT} ({\bf s}_i | \boldsymbol{\theta })\;. \end{equation} \tag{ 53 }$$

Let us define the "average" PDF as the product of the prior distribution of the parameters multiplied by the geometric mean of the likelihoods which describe these measurements:

$$\begin{equation} p_{\rm ave}(\boldsymbol{\theta } | {\bf s}) \equiv {\cal N} \, p^{(0)} {\cal L}^{\rm joint}_{\rm TOT}({\bf s}_1, {\bf s}_2, \ldots,\,\ {\bf s}_M| \boldsymbol{\theta })^{1/M}\;. \end{equation} \tag{ 54 }$$

Expanding this definition, we find

$$\begin{eqnarray} p_{\rm ave} (\boldsymbol{\theta } | \bf {s}) & \equiv & {\cal N} \, p^{(0)} \prod _{i=1}^{M}\, \left[ {\cal L}_{\rm TOT} ({\bf s}_i | \boldsymbol{\theta }) \right]^{1/M}\;, \end{eqnarray} \tag{ 55 }$$

where the subscript i denotes the ith noise realization in our set of M observations. The "ensemble average likelihood function" can in turn be expanded as

$$\begin{eqnarray} \prod _{i=1}^{M} [{\cal L}_{\rm TOT} ({\bf s}_i | {\boldsymbol{\theta }})]^{1/M} & = & \prod _{a} \prod _{i=1}^{M} [{\cal L}_{a} (s_{a,i} \, | \,{\boldsymbol{\theta }})]^{1/M}\nonumber\\ &=& \prod _{a} \prod _{i=1}^{M} e^{ - (h_a({\boldsymbol{\theta }}) - s_{a,i} \, | \, h_a({\boldsymbol{\theta }}) - s_{a,i})/2M } \nonumber \\ & = & \prod _{a} e^{-(h_a({\boldsymbol{\theta }}) - h_a({\hat{\boldsymbol\theta }}) \, | \, h_a({\boldsymbol{\theta }}) - h_a ({\hat{\boldsymbol\theta }}))/2} \nonumber \\ &&\times \prod _{i=1}^M \exp \left[\frac{1}{M}(n_{a,i} \, | h_{a}(\boldsymbol{\theta }) - h_{a}({\hat{\boldsymbol\theta }})) \right] \nonumber \\ &&\times \prod _{i=1}^M \exp \left[-\frac{1}{2M}(n_{a,i} \, | n_{a,i})\right].\quad \end{eqnarray} \tag{ 56 }$$

By taking M to be large, the last two lines of Equation (56) can be evaluated as follows:

$$\begin{eqnarray} &&\prod _{i=1}^M \exp \left[\frac{1}{M}(n_{a,i} \, | h_{a}(\boldsymbol{\theta }) - h_{a}({\hat{\boldsymbol\theta }})) \right]\nonumber\\ &&\qquad = \exp \left[\frac{1}{M}\sum _{i = 1}^M (n_{a,i} \, | h_{a}(\boldsymbol{\theta }) - h_{a}({\hat{\boldsymbol\theta }})) \right] \nonumber \\ &&\qquad\simeq \exp [\langle (n_{a} \, | h_{a}(\boldsymbol{\theta }) - h_{a}({\hat{\boldsymbol\theta }}))\rangle] = 1. \end{eqnarray} \tag{ 57 }$$

Here, 〈...〉 denotes an ensemble average over noise realizations (cf. Section 3.1), and we have used the fact that our noise has zero mean. Similarly, we find

$$\begin{eqnarray} \prod _{i=1}^M \exp \left[-\frac{1}{2M}(n_{a,i} \, | n_{a,i})\right] &=& \exp \left[-\frac{1}{2M}\sum _{i = 1}^M (n_{a,i} \, | n_{a,i}) \right] \nonumber \\ &\simeq & \exp \left[-\frac{1}{2}\langle (n_{a} \, | n_{a})\rangle \right] \nonumber \\ &=& e^{-1}\;. \end{eqnarray} \tag{ 58 }$$

This uses 〈(n_a|n_a)〉 = 2, which can be proved using the noise properties (35), (36), and (37).

Putting all this together, we finally find

$$\begin{equation} p_{\rm ave}(\boldsymbol{\theta } | {\bf s}) = {\cal N} p^{0}(\boldsymbol{\theta }) \prod _{a} e^{ - (h_a({\boldsymbol{\theta }}) - h_a ({\hat{\boldsymbol\theta }}) \, | \, h_a({\boldsymbol{\theta }}) - h_a ({\hat{\boldsymbol\theta }}))/2 }\;, \end{equation} \tag{ 59 }$$

where we have absorbed e⁻¹ into the normalization ${\cal N}$. The posterior PDF, averaged over noise realizations, is simply obtained by evaluating Equation (39) with the noise n set to zero.

Validation of our MCMC results requires comparing to work which goes beyond the Gaussian approximation and Fisher matrix estimators. In Section IVD of CF94, Cutler & Flanagan investigate effects that are nonlinear in 1/S/N. They show that such effects have a significant impact on distance measurement accuracies for low S/N. In particular, they find that Fisher-based estimates understate distance measurement errors for a network of two LIGO detectors and Virgo.

Because they go beyond a Fisher matrix analysis, the results of CF94 are useful for comparing with our results. Their paper is also useful in that they take source position to be known. Our approach is sufficiently different from CF94 that we do not expect perfect agreement, however. The most important difference is that we directly map out the posterior PDF and compute sample averages using Equations (48) and (49), for the full parameter set $\lbrace {\cal M}_z, \mu _z, D_L, \cos \iota, \psi, t_c, \Phi _c \rbrace$. In contrast, CF94 estimate measurement errors only for D_L, using an approximation on an analytic Bayesian derivation of the marginalized PDF for D_L. Specifically, Cutler & Flanagan expand the exponential factor in Equation (39) beyond second order in terms of some "best-fit" maximum-likelihood parameters. Their approximation treats strong correlations between the parameters D_L and cos ι that are nonlinear in 1/S/N. However, other correlations between D_L and (ψ, ϕ_c) are only considered to linear order. They obtain an analytic expression for the posterior PDF of the variables D_L and cos ι in terms of their "best-fit" maximum-likelihood values $\tilde{D}_L$ and $\cos \tilde{\iota }$ (see Equation (4.57) of CF94). The marginalized 1-D posterior PDFs for D_L are then computed by numerically integrating over cos ι. The one-dimensional marginalized PDF we compute in parameter θ_i is

$$\begin{equation} p_{\rm marg}(\theta _i | {\bf s}) = \int \dots \int p(\boldsymbol{\theta } | {\bf s}) d\theta _1 \dots d\theta _{i-1}\; d\theta _{i+1} \dots d\theta _N, \end{equation} \tag{ 60 }$$

where $p(\boldsymbol{\theta } | \bf {s})$ is the posterior PDF given by Equation (39) and N is the number of dimensions of our parameter set.

In addition to this rather significant difference in techniques, there are some minor differences which also affect our comparison.

• 1.  
  We use the restricted 2PN waveform; CF94 use the leading "Newtonian, quadrupole" waveform that we used for pedagogical purposes in Section 1.2. Since distance is encoded in the waveform's amplitude, we do not expect that our use of a higher-order phase function will have a large impact. However, to avoid any easily circumvented mismatch, we adopt the Newtonian-quadrupole waveform for these comparisons. This waveform does not depend on reduced mass μ, so for the purpose of this comparison only, our parameter space is reduced from 7 to 6 dimensions.
• 2.  
  We use the projected advanced sensitivity noise curve shown in Figure 1; CF94 use an analytical form (their Equation (2.1)¹¹) based on the best-guess for what advanced sensitivity would achieve at the time of their analysis. Compared to the most recent projected sensitivity, their curve underestimates the noise at middle frequencies (∼40–150 Hz) and overestimates it at high frequencies (≳200 Hz). We adopt their noise curve for this comparison. Because of these differences, CF94 rather seriously overestimates the S/N for NS–NS inspiral. Using their noise curve, the average S/N for the binary analyzed in their Figure 10 is 12.4¹²; using our up-to-date model for advanced LIGO, it is 5.8. As such, the reader should view the numbers in this section of our analysis as useful only for validation purposes.
• 3.  
  The two analyses use different priors. As extensively discussed in Section 3.3, we set uniform priors on the chirp mass ${\cal M}_z$, on the time t_c and phase Φ_c at coalescence, and on the polarization phase ψ. For this comparison, we assume isotropic emission and set a flat prior on cos ι. We assume that our sources are uniformly distributed in constant comoving volume. However, our detection threshold depends on the total network S/N, and effectively sets a joint prior on source inclination and distance. CF94 use a prior distribution only for the set {D_L, cos ι, ψ, Φ_c} that is flat in polarization phase, coalescence phase, and inclination. They assume a prior that is uniform in volume, but that cuts off the distribution at a distance D_L,max ≃ 6.5 Gpc.

Our goal here is to reproduce the 1-D marginalized posterior PDF in D_L for the binary shown in Figure 10 of CF94. We call this system the "CF binary." Each NS in the CF binary has m_z = 1.4 M_☉ and sky position (θ, ϕ) = (50°, 276°); the detector network comprises LIGO Hanford, LIGO Livingston, and Virgo. CF94 report the "best-fit" maximum-likelihood values ($\tilde{D}_L$, $\cos \tilde{\iota }$, $\tilde{\Psi }$) to be (432 Mpc, 0.31, 1015), where Ψ = ψ + Δψ(n), and where Δψ(n) depends on the preferred basis of e^× and e^× set by the detector network (see Equations (4.23)–(4.25) of CF94¹³). To compare our distribution with theirs, we assume that ${\hat{\boldsymbol\theta }} = {\tilde{\boldsymbol\theta }}_{\rm ML}$ for the purpose of computing the likelihood function ${\cal L}(\boldsymbol{\theta } | {\bf s})$. This is a reasonable assumption when the priors are uniform over the relevant parameter space. As already mentioned, for this comparison we use their advanced detector noise curve and the Newtonian-quadrupole waveform. Finally, we interpret the solid curve in Figure 10 of CF94 as the marginalized 1-D posterior PDF in D_L for an average of posterior PDF of parameters (given an ensemble of many noisy observations for a particular event). We compute the average PDF as described in Section 3.4, and then marginalize over all parameters except D_L, using Equation (60).

The left-hand panels of Figure 5 show the resulting one-dimensional marginalized PDF in D_L and cos ι. Its shape has a broad structure not dissimilar to the solid curve shown in Figure 10 of CF94: The distribution has a small bump near D_L ≈ 460 Mpc, a main peak at D_L ≈ 700 Mpc, and extends out to roughly 1 Gpc. Because of the broad shape, the Bayes mean ($\tilde{D}_{L,\rm {BAYES}} = 694\,{\rm Mpc}$) is significantly different from both the true value ($\hat{D}_L = 432\,{\rm Mpc}$ in our calculation) and from the maximum likelihood ($\tilde{D}_{L,\rm {ML}} = 495\,{\rm Mpc}$). Thanks to the marginalization, the peak of this curve does not coincide with the maximum likelihood.

Zoom In Zoom Out Reset image size

We further determine the 2-D marginalized posterior PDFs in D_L and cos ι for the CF binary. Figure 5 illustrates directly the very strong degeneracy between these parameters, as expected from the form of Equations (7) and (8), as well as from earlier works (e.g., Marković 1993, CF94). It is worth noting that, as CF94 comment, this binary is measured particularly poorly. This is largely due to the fact that one polarization is measured far better than the other, so that the D_L–cos ι degeneracy is essentially unbroken. This degeneracy is responsible for the characteristic tail to large D_L we find in the 1-D marginalized posterior PDF in D_L, $p(D_L | \bf {s})$, which we investigate further in the following section.

We now examine how well we measure D_L as a function of distance to the CF binary and the properties of the GW detector network. Figures 6 and 7 show the 1-D and 2-D marginalized posterior PDFs in D_L and cos ι for the CF binary at $\hat{D}_L = \lbrace 100$, 200, 300, 400, 500, 600} Mpc. For all these cases, we keep the binary's sky position, inclination, and polarization angle fixed as in Section 4.1. The average network S/Ns we find for these six cases are (going from $\hat{D}_L = 100\,{\rm Mpc}$ to 600 Mpc) 53.6, 26.8, 17.9, 13.4, 10.7, and 8.9 (scaling as $1/\hat{D}_L$). Interestingly, the marginalized PDFs for both distance and cos ι shown in Figures 6 and 7 have fairly Gaussian shapes for $\hat{D}_L = 100$ and 200 Mpc, but have very non-Gaussian shapes for $\hat{D}_L \ge 300\,{\rm Mpc}$. This can be considered "anecdotal" evidence that the Gaussian approximation for the posterior PDF breaks down at S/N ≲ 25 or so, at least for this case. For lower S/N, the degeneracy between cos ι and D_L becomes so severe that the 1-D errors on these parameters become quite large.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Next, we examine measurement accuracy versus detector network. For the CF binary, adding detectors does not substantially increase the total S/N. We increase the average total S/N from 12.4 to 14.6 (adding only AIGO), to 12.4 (adding only LCGT; its contribution is so small that the change is insignificant to the stated precision), or to 14.7 (adding both AIGO and LCGT). The average S/N in our detectors is 8.23 for LIGO-Hanford, 8.84 for LIGO-Livingston, 2.91 for Virgo, 8.71 for AIGO, and 1.1 for LCGT. This pathology is an example of a fairly general trend that we see; it is common for the S/N to be quite low in one or more detectors.

In the case of the CF binary, we find that adding detectors does not improve the measurement enough to break the D_L–cos ι degeneracy. The marginalized PDFs as functions of D_L and cos ι remain very similar to Figure 5, so we do not show them. As a consequence, even with additional detectors, the distance errors remain large and biased. The bias is because we tend to find cos ι to be larger than the true (relatively edge-on) value (cf. the lower left-hand panel of Figure 5). Thanks to the D_L–cos ι degeneracy, we likewise overestimate distance.

One of the primary results from the CF binary analysis is a strong degeneracy between cos ι and D_L. As Figure 5 shows, this results in a tail to large distance in the one-dimensional marginalized posterior PDF $p(D_L | \bf {s})$, with a Bayes mean $\tilde{D}_L = 694\,{\rm Mpc}$ (compared to $\hat{D}_L = 432\,{\rm Mpc}$). Such a bias is of great concern for using binary sources as standard sirens.

The CF binary has $\cos \hat{\iota }= 0.31$, meaning that it is nearly edge-on to the line of sight. Hypothesizing that the large tails may be due to its nearly edge-on nature, we consider a complementary binary that is nearly face on: we fix all of the parameters to those used for the CF binary, except for the inclination, which we take to be $\cos \hat{\iota }= 0.98$. We call this test case the "face-on" CF binary. Changing to a more nearly face-on situation substantially augments the measured S/N; the average S/N for the face-on CF binary measured by the LIGO/Virgo base network is 24.3 (versus 12.4 for the CF binary). We thus expect some improvement simply owing to the stronger signal.

Figure 8 shows the 1-D and 2-D marginalized posterior PDFs in D_L and cos ι. As expected, these distributions are complementary to those we found for the CF binary. In particular, the peak of the 1-D marginalized posterior PDF in D_L is shifted to lower values in D_L, and the Bayes mean is much closer to the true value: $\tilde{D}_L = 376.3\,{\rm Mpc}$. The shape of the 1-D marginalized posterior PDF in cos ι is abruptly cut off by the upper bound of the physical prior cos ι ⩽ 1, and the tail extends to lower distances (the opposite of the CF binary). The Bayes mean for the inclination is $\cos \tilde{\iota }= 0.83$.

Zoom In Zoom Out Reset image size

Just as we varied distance and detector network for the CF binary, we also do so for the face-on CF binary, with very similar results. In particular, varying network has little impact on the marginalized 1-D PDFs in D_L and cos ι. Varying distance, we find that the marginalized 1-D PDFs are nearly Gaussian in shape for small distances, but become significantly skewed (similar to the left-hand panels of Figure 8) when $\hat{D}_L > 200$ Mpc. The distributions in cos ι are particularly skewed thanks to the hard cutoff at cos ι = 1. Interestingly, in this case we tend to infer a value of cos ι that is smaller than the true value. We likewise find a Bayes mean $\tilde{D}_L$ that is smaller than $\hat{D}_L$.

The main result from our testing is that the posterior PDFs we find have rather long tails, with strong correlations between cos ι and D_L. Except for cases with very high S/N, the one-dimensional marginalized posterior PDF in cos ι tends to be rather broad. The Bayes mean for cos ι thus typically suggests that a binary is at intermediate inclination. As such, we tend to underestimate cos ι for nearly face-on binaries, and to overestimate it for nearly edge-on binaries. Overcoming this limitation requires us to either break the D_L–cos ι degeneracy (such as by setting a prior on binary inclination), or by measuring a population of coalescences. Measuring a population will make it possible to sample a wide range of the cos ι distribution, so that the event-by-event bias is averaged away in the sample.

We begin by imagining a population of 600 detected NS–NS binaries, either isotropically distributed in inclination angle or from our beamed subsample, using a network with all five detectors. Figure 9 shows scatter plots of the distance measurement accuracies for our unbeamed (blue crosses) and beamed events (black dots), with each panel corresponding to a different detector network. The distance measurement error is defined as the ratio of the rms measurement error with the true value¹⁴$\hat{D}_L$:

$$\begin{equation} \frac{\Delta D_L}{\hat{D}_L} = \frac{\sqrt{\Sigma ^{D_L D_L}}}{\hat{D}_{L}}\;. \end{equation} \tag{ 62 }$$

$\Sigma ^{D_L D_L}$ is computed using Equation (51). We emphasize some general trends in Figure 9 which are particularly relevant to standard sirens.

• 1.  
  The unbeamed total sample and the beamed subsample separate into two distinct distributions. As anticipated, the beamed subsample improves measurement errors in D_L significantly, by greater than a factor of two or more. This is due to the beaming prior, which constrains the inclination angle, cos ι, to ∼3%, thereby breaking the strong D_L–cos ι degeneracy. By contrast, when no beaming prior is assumed, we find absolute errors of 0.1–0.3 in cos ι for the majority of events. The strong D_L–cos ι degeneracy then increases the distance errors. A significant fraction of binaries randomly selected from our sample have $0.5 \lesssim |\cos \hat{\iota }| < 1$. As discussed in Section 3.2, this is due to the S/N selection criterion: at fixed distance, face-on binaries are louder and tend to be preferred.
• 2.  
  Beamed subsample scalings. We fit linear scalings to our beamed subsample:$\Delta D_L/\hat{D}_L \simeq \hat{D}_L/(2.15\,\rm {Gpc})$ for LIGO + Virgo$\Delta D_L/\hat{D}_L \simeq \hat{D}_L/(2.71\,\rm {Gpc})$ for LIGO + Virgo + AIGO$\Delta D_L/\hat{D}_L \simeq \hat{D}_L/(2.38\,\rm {Gpc})$ for LIGO + Virgo + LCGT$\Delta D_L/\hat{D}_L \simeq \hat{D}_L/(2.82\,\rm {Gpc})$ for LIGO + Virgo + AIGO + LCGT.
• 3.  
  When isotropic emission is assumed, we find a large scatter in distance measurement errors for all events, irrespective of network and true distance. We find much less scatter when we assume a beaming prior. This is illustrated very clearly by the upper-right panel of Figure 9. In that panel, we show the scatter of distance measurement error versus true distance for the LIGO, Virgo, AIGO detector network, comparing to the Fisher-matrix-derived linear scaling trend found in DHHJ06. For the unbeamed case, our current results scatter around the linear trend; for the beamed case, most events lie fairly close to the trend. This demonstrates starkly the failure of Fisher methods to estimate distance accuracy, especially when we cannot set a beaming prior.
• 4.  
  Adding detectors to the network considerably increases the number of detected binaries, but does not significantly improve the accuracy with which those binaries are measured. The increase we see in the number of detected binaries is particularly significant for GW-SHB standard sirens. For instance, an important application is mapping out the posterior PDF for the Hubble constant, H₀. As the number of events increases, the resulting joint posterior PDF in H₀ will become increasingly well constrained. Additional detectors also increase the distance to which binaries can be detected. This can be seen in Figure 9: for the LIGO and Virgo network, our detected events extend to $\hat{D}_L \sim 600\,{\rm Mpc}$; the larger networks all go somewhat beyond this. Interestingly, networks which include the AIGO detector seem to reach somewhat farther out.

Zoom In Zoom Out Reset image size

It is perhaps disappointing that increasing the number of detectors does not greatly improve measurement accuracy. We believe this is due to two effects. First, a larger network tends to detect more weak signals. These additional binaries are poorly constrained. Second, the principal limitation to distance measurement is the D_L–cos ι degeneracy. A substantial improvement in distance accuracy on individual events would require breaking this degeneracy. We find that adding detectors does not do this, but the beaming prior does.

We now repeat the preceding analysis for 600 detected NS–BH binaries. Figure 10 shows scatter plots of measurement accuracies for unbeamed and beamed NS–BH binaries. We find similar trends in the NS–NS case.

• 1.  
  The unbeamed and beamed samples separate into two distinct distributions. Note, however, that outliers exist in measurement errors at high D_L for several beamed events for all networks. This is not too surprising, given that we expect beamed sources at higher luminosity distances and lower S/N. Such events are more likely to deviate from the linear relationship predicted by the Fisher matrix.
• 2.  
  We see substantial scatter in distance measurement, particularly when isotropic emission is assumed. As with the NS–NS case, the scatter is not as severe when we assume beaming, and in that case lies fairly close to a linear trend, as would be predicted by a Fisher matrix. This trend is shallower in slope than for NS–NS binaries, thanks to the larger mass of the system.
• 3.  
  We do not see substantial improvement in distance measurement as we increase the detector network. As with NS–NS binaries, adding detectors increases the range of the network; AIGO appears to particularly add events at large $\hat{D}_L$ (for both the isotropic and beamed samples). However, adding detectors does not break the fundamental D_L–cos ι degeneracy, and does not improve errors. From our full posterior PDFs, we find absolute errors of 0.1–0.3 in cos ι, which is very similar to the NS–NS case.
• 4.  
  Beamed subsample scalings. The linear scalings for our beamed subsample are$\Delta D_L/\hat{D}_L \simeq \hat{D}_L/(4.83\,\rm {Gpc})$ for LIGO + Virgo$\Delta D_L/\hat{D}_L \simeq \hat{D}_L/(6.14\,\rm {Gpc})$ for LIGO + Virgo + AIGO$\Delta D_L/\hat{D}_L \simeq \hat{D}_L/(5.20\,\rm {Gpc})$ for LIGO + Virgo + LCGT$\Delta D_L/\hat{D}_L \simeq \hat{D}_L/(6.76\,\rm {Gpc})$ for LIGO + Virgo + AIGO + LCGT.

Zoom In Zoom Out Reset image size