The NINJA-2 catalog of hybrid post-Newtonian/numerical-relativity waveforms for non-precessing black-hole binaries

A new generation of laser interferometric gravitational-wave detectors (Advanced LIGO (aLIGO) [1–3], Advanced Virgo [4, 5] and LCGT [6]) is currently under construction. These second-generation detectors will have an order of magnitude increase in sensitivity over first-generation instruments and will be sensitive to a broader range of gravitational-wave frequencies. One of the most widely anticipated sources for this global network of observatories is the inspiral, merger and ringdown of a binary containing two black holes [7]. Detection of such a binary black-hole coalescence will allow astronomers and astrophysicists to directly observe the physics of black-hole spacetimes and to explore the strong-field conditions of Einstein's theory of general relativity [8].

The ability of gravitational-wave astronomers to use the new generation of observatories to detect and study binary black-hole coalescence depends on the quality of search and source-parameter measurement algorithms. These algorithms rely on the physical accuracy of the underlying theoretical waveform models. Developing and testing the algorithms required to achieve the goals of gravitational-wave astronomy demands close interaction between the source-modeling and data-analysis communities. The numerical injection analysis (NINJA) project was created in 2008 to bring these communities together and to use the recent advances in numerical relativity (NR) [9] to test analysis pipelines by adding physically realistic signals to detector noise in software. We describe such additions of signals into noise as 'injections'.

The first NINJA project (NINJA-1) [10] considered a total of 23 numerical waveforms, which were injected into Gaussian noise colored with the frequency sensitivity of first-generation detectors. These data were analyzed by nine data-analysis groups using both search and parameter-estimation algorithms. However, there were two major limitations to the NINJA-1 analysis. First, to encourage broad participation, no length or accuracy requirements were placed on the numerical waveforms. Consequently, many of these waveforms were too short to inject over an astrophysically interesting mass range without introducing artifacts into the data. The lowest mass binary considered in NINJA-1 had a total mass of 35 M_⊙, whereas the mass of black holes could be below 5 M_⊙ [11, 12]. In NINJA-1, the waveforms were only inspected for obvious, pathological errors and no cross-checks were performed between the submitted waveforms, and therefore, it was difficult to assess the physical fidelity of the results. Second, the NINJA-1 data set contained stationary noise with the simulated signals already injected into the data. The data set contained only 126 simulated signals, which precluded detailed statistical studies of the effectiveness of search and parameter-estimation algorithms. Finally, since the data set lacked the non-Gaussian noise transients present in real detector data, it was not possible to fully explore the response of the algorithms in a real search scenario. Despite these limitations, NINJA-1 successfully removed a number of barriers to collaboration between the source-modeling and data-analysis communities and demonstrated where further work is needed. The goal of the second NINJA project (NINJA-2) is to address the deficiencies of NINJA-1 and to perform a systematic test of the efficacy of data-analysis pipelines in real-world situations in preparation for aLIGO and Virgo.

This paper reports on the improvements we have made to the NINJA analysis to address the first of the limitations described above—the accuracy of the numerical waveforms. We present the NINJA-2 waveform catalog and describe the results of the procedures we have used to validate these data. NINJA-2 places requirements on the accuracy and length of the contributed waveforms and we have performed systematic cross-checks of the submitted waveforms. Each binary black-hole simulation in the NINJA-2 catalog must include at least five orbits of usable data before merger, i.e. neglecting the initial burst of junk radiation. The NR waveform amplitude should be accurate to within 5%, and the phase (as a function of gravitational-wave frequency) should have an accumulated uncertainty over the entire inspiral, merger and ringdown (of the numerical simulation) of no more than 0.5 rad. We also require that numerical simulations are 'hybridized' to post-Newtonian (pN) waveforms so that the resulting waveforms contain enough cycles to allow injections at M ⩾ 10 M_⊙ in early aLIGO data. The continued advances in numerical simulations have also allowed us to study a somewhat larger region of the signal parameter space; however, we have restricted our attention to non-precessing binaries for reasons discussed in section 2.

A subsequent paper will describe the results of using the NINJA-2 waveforms to study search and parameter-estimation algorithms in real detector data. Since data from the second-generation detectors are not yet available, the NINJA-2 analysis will use data from the first-generation detectors re-colored so that it has the frequency response expected in the first observing runs of the aLIGO detectors—referred to as 'early aLIGO' [13]. Similar noise curves will be used to simulate the Advanced Virgo detector. NINJA-2 analyses will use these noise models to ensure that existing algorithms are optimal when second-generation detectors come online in ∼2015. The results in this paper use the early aLIGO sensitivity curve (cf figure 1) to study the accuracy of the submitted waveforms. The ultimate sensitivity of aLIGO is expected to be significantly better than this curve. To allow the waveforms to be used in studies using more sensitive noise curves, we have also performed accuracy studies using the aLIGO zero-detuned, high-power [14] sensitivity curve. Figure 1 shows the two aLIGO sensitivity curves, characterized by their amplitude spectral densities (ASD) overlaid with one of the contributed NINJA-2 waveforms. This figure demonstrates that hybridization is necessary to allow scaling of the numerical waveforms to astrophysically interesting masses, and a portion of this paper studies the hybridization methods used to construct the NINJA-2 waveforms.

Zoom In Zoom Out Reset image size

The remainder of this paper is organized as follows. Section 2 describes in more detail the accuracy requirements that we have placed on NINJA-2 simulations and presents an overview of the waveform catalog showing the regions of the binary black-hole parameter space covered. Section 3 gives an overview of the numerical methods used to construct the NR waveforms and section 4 describes the methods that we have used to hybridize the numerical simulations to pN waveforms. The pN waveforms themselves are summarized in the appendix. Section 5 describes the methods and results of the comparisons we have performed between the waveforms. Based on these comparisons, we judge the hybrid waveforms suitable for the NINJA-2 project. Section 6 summarizes our findings and suggests directions for future improvements of the catalog, in particular, the study of higher order modes in the waveforms.

Binary black holes formed from the evolution of massive stars are expected to have circularized before their gravitational-wave frequency reaches the sensitive band of ground-based detectors, and so we only consider circular (non-eccentric) binaries. In the NINJA-2 project, we further restrict our attention to binaries that do not undergo precession, i.e. where the spins of the black holes either vanish or are parallel (or anti-parallel) to the binary's orbital angular momentum. We do this for two reasons: (i) in trying to understand the complex phenomenology of the binary parameter space, we prefer to tackle first a simpler subset, which nonetheless captures the main features of binary inspiral and merger and (ii) the precessing-binary parameter space has been sampled by only a handful of numerical simulations. The NR community is currently exploring the space of precessing binaries, e.g., through the numerical relativity–analytical relativity (NR–AR) project [15]. Such waveforms will be used in future NINJA projects that explore precession.

The parameters of the black-hole binaries we consider in NINJA-2 are the mass of each black hole, m₁ and m₂, or equivalently, the total mass M = m₁ + m₂ and the mass ratio q = m₁/m₂, and the dimensionless spin magnitude of each black hole, χ₁ ≡ S₁/m²₁ and χ₂ = S₂/m²₂. The total mass sets the overall scale of the system and can be factored out to leave a three-dimensional parameter space of waveforms, {q, χ₁, χ₂}. (In data-analysis applications, the mass is factored back in, and the space of intrinsic parameters for non-precessing binaries is four dimensional.) Figure 2 shows the coverage of parameter space (details in section 3). The sampling is coarse; while the waveforms will provide invaluable information within the NINJA-2 project, we expect that ultimately a more uniform coverage of parameter space by a much larger number of configurations will be necessary.

Zoom In Zoom Out Reset image size

The NINJA-2 requirement of five pre-merger orbits is at the low end of estimates of sufficient waveform lengths for the construction of accurate hybrid pN–NR waveforms, as discussed in [16–20], but we expect these to be acceptable for the goals of the NINJA-2 project. The 5% amplitude and 0.5 rad phase accuracy requirements were formulated with typical current waveforms in mind, e.g., those studied in the Samurai project [21] and studies performed in preparation for the NR–AR collaboration project. These requirements are consistent for waveforms of similar lengths but may not be directly applicable to much longer waveforms. For example, in the 25-orbit SpEC simulations with dimensionless spins χ_i = 0.97, the highest- and second-highest-resolution data differ by roughly 0.6 rad at merger [22]. Note that because this phase-error accumulates over 20 additional inspiral orbits, this waveform would easily satisfy the NINJA-2 phase requirement if it were truncated to minimally meet the NINJA-2 length requirement (although such a truncation would decrease the accuracy of the hybridized waveform).

We require that the hybridized waveforms in the NINJA-2 catalog are long enough to span the sensitivity bands of the aLIGO and Virgo detectors in their early operation. Specifically, when rescaled to 10 M_⊙, the hybrids must begin with a gravitational wave frequency of 20 Hz or lower, i.e. a starting GW frequency of Mω ⩽ 0.006. This requires extending the NR waveforms to lower frequencies (i.e. more inspiral cycles) by attaching a pN inspiral waveform onto the early portion of the NR waveform to produce a hybrid pN–NR waveform. We require that the hybridization be performed at a gravitational-wave frequency of Mω₂₂ ⩽ 0.075, where Mω₂₂ is the frequency of the (ℓ, m) = (2, ±2) harmonic. This is a pragmatic choice, consistent with our requirement of five pre-merger orbits and supported by studies of the physical accuracy of the late-inspiral phase evolution of pN waveforms for a selection of points in the non-precessing-binary parameter space [23–27]. In practice, hybridization fits were performed over a frequency range as summarized in section 4, and the average frequency, and frequency of the average time of the fitting interval were always chosen below Mω₂₂ ⩽ 0.075, with two exceptions as given in table 1: the non-spinning Llama waveforms at mass ratios q = 1, 2. As shown in figure 7, these do however show excellent overlaps with comparison waveforms.

The waveforms were submitted with the complex GW strain function h₊ − ih_× decomposed into modes using spin-weighted spherical harmonics ⁻²Y_ℓm of weight s = −2. Although most of the power is in the (ℓ, m) = (2, ±2) modes, we encouraged (but did not require) the submission of additional subdominant modes. The accuracy studies in this paper focus on the (ℓ, m) = (2, 2) mode; further work is required to study the accuracy of the contributed subdominant modes. A total of 63 waveforms from 8 groups were contributed to the NINJA-2 catalog. There are 40 distinct numerical waveforms; some of these waveforms have been hybridized with multiple pN waveforms. The NINJA-2 catalog is summarized in table 1, and a map of the parameter values is shown in figure 2. In the next section, we describe in more detail the numerical methods used to generate these waveforms and present additional plots in figures 3 and 4.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

3.1. Summary of contributions

3.2. Data format of contributions

3.3. Initial data and eccentricity

3.4. Gravitational-wave extraction

3.5. Numerical methods, accuracy and sources of error

While pN methods accurately approximate GW signals throughout early inspiral, they become increasingly unreliable toward late inspiral and merger [122]. NR simulations are capable of accurately computing inspiral, merger and ringdown portions of binary coalescence [59, 78, 97, 123], but are computationally too expensive to extend far into the inspiral regime [32]. Hybrid waveforms are the result of smoothly blending together pN and NR waveforms to form a waveform that covers the full BBH dynamics. For the NINJA-2 project, each NR group has produced their own hybrid waveform after ensuring that the pN portions all agree. We expect some systematic errors resulting from errors in the pN approximation, the choice of blending region and the hybridization method [16, 18–20, 32, 124, 125]. The NINJA-2 hybrid waveforms do not contain EOB extensions of pN approximants, although these have also been used to model complete waveforms (e.g., [19, 126, 127]).

Hybridization uses least-squares fits to determine the extrinsic parameters for the pN waveform [125, 128, 129]. In general, this is accomplished by evaluating

$$\begin{eqnarray} &&\delta ( \vec{u}, a ) = \min _{\lbrace \vec{u},a\rbrace } \int ^{\xi _{1}}_{\xi _{2}}{|\Upsilon _{{\rm pN}}(\xi ,\vec{u}) - a\Upsilon _{{\rm NR}} (\xi ,\vec{u}_{0})|^2\, \mathrm{d}\xi }, \end{eqnarray} \tag{ 1 }$$

where ϒ represents waveform data relating to strain [e.g., h(t) = h₊(t) − i h_×(t), $\arg [h(t)]$ or $\tilde{h}(f)$]. If ϒ is derived from the time domain, then ξ = t; if ϒ is in the frequency domain, then ξ = f. For either case, [ξ₁, ξ₂], chosen within the domain of both the pN and NR data sets, defines the integration interval and, in most cases, the blending region. The vector $\vec{u}$ denotes the set of pN parameters over which the fitting is performed. For example, $\vec{u} = (t_{{\rm shift}},\phi _{{\rm shift}},q)$ corresponds to adjusting time and phase shift and the mass ratio of the pN waveform to match the NR waveform. The best-fit parameters are denoted by $\vec{u}^\ast$. The amplitude scaling factor a is often fixed to a = 1, but may be included in the fitting parameters [128]. Finally, in the limit ξ₁ → ξ₂, this procedure reduces to enforcing equality of ϒ_pM and ϒ_NR at ξ₁ = ξ₂, as well as the equality of the first derivative.

More explicitly, hybridization may be performed via the following algorithm.

• (i)  
  Choose [ξ₁, ξ₂] within the pN and NR data sets. Ideally, [ξ₁, ξ₂] is sufficiently early so that both pN and NR sets should be accurate.
• (ii)  
  Evaluate equation (1); apply $\lbrace \vec{u}^\ast , a^\ast \rbrace$ to the pN data set, resulting in ϒ*_pN. Measure error quantities relating to fit.
• (iii)  
  If desired, adjust [ξ₁, ξ₂] and iterate (i) and (ii), to find a preferred interval [ξ*₁, ξ₂*].
• (iv)  
  Defining a monotonic function z(ξ), such that z(ξ < ξ_a) = 0 and z(ξ > ξ_b) = 1, the hybrid is given by

  $$\begin{equation} \Upsilon _{{\rm Hyb}}(\xi ) = [1-z(\xi )]\, \Upsilon _{{\rm pN}}^\ast + a^{\ast }\, z(\xi )\, \Upsilon _{{\rm NR}} . \end{equation} \tag{ 2 }$$

  Note that the transition region [ξ_a, ξ_b] is generally taken to be a sub-interval of [ξ*₁, ξ₂*], sometimes consisting of a single point.

This general formalism leaves many decisions open, and the following specific choices were made during hybridization.

• For the SpEC waveforms, the integrand of equation (1) was taken to be the square of the phase difference only [129], $\Upsilon = \arg [h(t)]$, and maximation was performed over t_shift and ϕ_shift only, without any adjustments to the amplitude (a* ≡ 1). The time-interval [t₁, t₂] was chosen to correspond to the frequencies listed in table 1, without any iterative adjustments of t₁ and t₂.
• The RIT waveforms similarly use $\Upsilon = \arg [h(t)]$ for equation (1), but employ the limit as t₁ → t₂ to determine $\vec{u} = (t_{{\rm shift}},\phi _{{\rm shift}})$ at Mω = 0.075. Then, the transition function is z_RIT(t) = x(t)³ [6x(t)² − 15x(t) + 10], where x(t) := (t − t₁)/(t₂ − t₁) for a finite interval (t₁ ≠ t₂), guaranteeing C² behavior at t = t₁ and t₂ [130].
• For the Lean waveforms hybridization is performed in a similar way [131], using the transition function z(t) = 70x⁹ − 315x⁸ + 540x⁷ − 420x⁶ + 126x⁵and individual mode amplitudes of the pN waveforms are rescaled such that their average over the matching window agrees with the numerical result.
• The UIUC waveforms use ϒ = h(t), a = 1, using as free parameters the initial pN phase and orbital angular frequency. Equation (2) is then used to construct the final hybrid, with z(t) = (t − t₁)/(t₂ − t₁) as in [132]. The time interval [t₁, t₂] was chosen to correspond to the UIUC frequencies listed in table 1, without any iterative adjustments of t₁ and t₂.
• The GATech hybridization follows [129] and is done in the time domain with ϒ = h(t) and $\vec{u} = (t_{{\rm shift}}, \phi _{{\rm shift}})$. Equation (1) is evaluated over $\lbrace \vec{u},a\rbrace ,$ and then, equation (2) is used to construct the final hybrid, with z(t) = (t − t₁)/(t₂ − t₁). The fitting intervals are given in table 1.

Figures 3 and 4 show exemplary plots of the resulting hybrid waveforms. For aligned spins, orbital hang-up extends the inspiral to smaller separation and higher frequencies. This is apparent in figure 3 that the last ten GW cycles (as indicated by the small circles) take less time and in figure 4 by the shift toward higher frequencies.

Each NR group verified that their waveforms met the minimum NINJA-2 requirements before submission, as described in section 1. The minimum-five-orbit requirement was easily verified by inspection, and the amplitude and phase uncertainties were estimated by convergence tests with respect to numerical resolution and waveform-extraction radius. Once submitted, a series of checks were performed in order to validate the waveforms against each other. In the first stage, the pN expressions and codes were compared against each other and the literature. This resulted in a set of codes in various languages producing waveforms that agree well in both phase and amplitude (see the appendix).

5.1. Time-domain and frequency-domain checks

5.2. Overlap comparisons

In this check, the waveforms were compared against each other using matched-filtering techniques. The inner product between two real waveforms s₁(t) and s₂(t) is defined as

$$\begin{equation} \left({s_1}|{s_2}\right) = 4\, \Re \int _{0}^\infty {\rm d}f\, \frac{\tilde{s}_1(f) \tilde{s}_2^\star (f)}{S_n(f)}, \end{equation} \tag{ 3 }$$

where S_n(f) is the power spectral density (see figure 1). The overlap is then obtained by normalization and maximization over relative time and phase shifts, Δt and Δϕ:

$$\begin{equation} \left\langle {s_1}|{s_2}\right\rangle :=\max _{\Delta t, \Delta \phi } \frac{ \left({s_1}|{s_2}\right) }{ \sqrt{ \left({s_1}|{s_1}\right) \left({s_2}|{s_2}\right) }}. \end{equation} \tag{ 4 }$$

There is an important subtlety involved in calculating these overlaps. Usually, all possible time shifts are evaluated simultaneously by a suitable combination of fast Fourier transforms. This results in overlaps at time shifts at discrete values of Δt spaced by the inverse sampling frequency. The result of time maximization is then taken to be the maximum of this discrete time series. When comparing two very similar waveforms, such as two NR simulations of the same system, the overlap function becomes very sharply peaked. In units, where G = c = 1, 1M_⊙ = 4.93 × 10⁻⁶ s, and in these units, time shifts of well below 1M can lead to significant changes in the overlap. It is therefore imperative that the sample rate used in calculating the overlap be large enough to find the true maximum.

This issue is demonstrated in figure 5. The left panel shows the overlap function resulting from the comparison of two equal-mass, non-spinning waveforms sampled at four different rates. At 4096 Hz, the peak of the function is missed and the overlap is underestimated. The consequence of this is illustrated in the right panel of figure 5, which shows that the overlap as a function of mass exhibits oscillations. As the mass changes, the phase shift Δϕ that maximizes equation (3) also changes. For some masses, the optimal phase shift will occur at one of the discretely sampled time shifts, and the results will be correct. For other masses, the true extremum will occur between discretely sampled time shifts (as for 4096 Hz in the left panel), and the computed overlap will be erroneously too low. Henceforth, in this paper, all overlaps are calculated at 32 768 Hz. The LIGO/Virgo matched filter searches operate on data at 4096 Hz; however, this issue is not a problem in these searches for reasons that can be seen in these two images. The loss of overlap by undersampling is no larger than 0.2%. The search utilizes a bank of templates that discretizes the parameter space. In constructing the template bank, we have already incurred a potential loss of SNR of 3%, and this effect is therefore negligible.

Zoom In Zoom Out Reset image size

Before discussing overlaps between the submitted NINJA hybrid waveforms, we first present comparisons between time-domain pN waveforms of the kind used to construct hybrid waveforms. (See the appendix for a summary of the pN approximants used.) These waveforms terminate at Mω = 0.136, 0.114, 0.069, and 0.135 for T1, T2, T3 and T4, respectively. At M = 10 M_⊙, these correspond to termination frequencies of 435, 369, 222 and 439 Hz, respectively. The results of our overlap calculations are shown in figure 6. In [133], the authors made a far more extensive comparison between pN approximants, although that analysis differs in several important aspects from what has been done here; notably, the overlaps are maximized over the intrinsic mass parameters of one waveform. However, the qualitative conclusion that T1 matches better with T2 and T4 than with T3 is consistent with our results. These results should be borne in mind when evaluating overlaps between hybrid waveforms using different pN approximants. In particular, overlaps at lower masses will be dominated by the influence of the two pN waverforms, and therefore, hybrid waveforms using T3 will have relatively low matches against hybrid waveforms using T1. Note also that T3 and T4 diverge from T1 as mass increases, this is precisely why we need to transition to NR waveforms.

Zoom In Zoom Out Reset image size

Let us now discuss the main results of this section, the overlap between different submitted hybrid waveforms. For six of the black-hole configurations listed in table 1, hybrids have been constructed independently for different numerical waveforms, and overlaps were computed between each pair of waveforms in a group. Below we are showing overlap plots for all six cases: the q = {1, 2} non-spinning waveform submissions are shown in figure 7, q = {3, 4} non-spinning waveforms in figure 8, and q = 1, χ_i = 0.4, 0.85 spinning waveforms in figure 9. We will discuss the q = {1, 2} non-spinning cases in more detail.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

At the high-mass end, where the NR portion of the waveform dominates the overlap, the overlaps approach 1. Since all these waveforms model the same physical system, this is the expected behavior. High overlaps at high masses indicate good agreement between different NR codes, and the results here are consistent with the detailed study of the q = 1 case that was performed in the Samurai project [21]. At the low-mass end, where the overlap is dominated by the pN portion of the waveform, the behavior is qualitatively as expected from figure 6. In particular, the overlap between the BAM hybrid using T1 and the SpEC hybrid using T4 is lower than the matches between the other hybrids, all of which use T4. In the region between ∼15–30 M_⊙, there are drops in the matches between hybrids using T4. In this mass range, the hybridization regions are passing through the most sensitive frequency band, and the mismatches are due in part to different choices of hybridization methods and parameters. The question of how many NR cycles are needed in order to produce a robust hybrid waveform is an area of active research [16–20].

If the approximants are the same, then the mismatches will depend only on the differences in (1) the hybridization methods, (2) the hybridization frequencies (and windows) and (3) the NR data. We have performed tests to verify that the mismatches due to these effects are very small. We have made overlap calculations using a white-noise PSD, integrated between 10 and 100 Hz, so the hybridization region can pass fully into and out of band as the mass varied between 10 and 100M_⊙ and found that the maximum mismatch due to hybridization is 0.05% for the GATech and SpEC equal-mass non-spinning hybrids. This suggests that the contribution of the hybridization to the mismatch is very small, which is consistent with the results in [16]. Mismatches for masses M ≳ 150M_⊙ will be due only to differences between the numerical data, and we find these to be 0.1%, which is consistent with the results for equal-mass non-spinning binaries in [21] and for q = 2 non-spinning binaries in [125].

The overlap plots discussed thus far do not yet address the accuracy required for detection of gravitational waves. Broadly, in order to claim a detection, a signal must match well against at least one template waveform used in a search. Several methods of assessing detection accuracy have been proposed [134, 135]; however, here we take the simple aproach that the waveforms are sufficient for further detection studies if the overlap is above the standard 0.97 threshold for a loss of no more than 10% of signals in a search. In some cases above, the overlaps are not above this threshold, but the appropriate quantity to evaluate to address this question is the overlap maximized over all of the physical parameters in a search. We now extend these overlap studies by maximizing over one physical parameter, the mass of one of the waveforms, as well as the time and phase. If the overlaps are now all over 0.97 (which they are), then they are acceptable for use in search-related studies.

These extended overlaps also provide insight into the 'parameter-estimation question', as the error in parameter recovery between two hybrid waveforms gives a rough lower bound on the errors that we may expect in recovering parameters of hybrid injections with search templates. In practice, parameter recovery is likely to be worse than this, as it will involve maximizing over several parameters in addition to total mass, in addition to differences between the waveforms. Example plots using the equal-mass, non-spinning MayaKranc waveform as the signal and BAM hybridized with two different pN approximants as the template are shown in figure 10. Note that in the right panel of figure 10 the minimum overlap is below 0.97. This minimum match will be even worse for q > 1 binaries [16], but if maximization is done over the other physical parameters (mass ratio and spin), then the overlap will increase to well above 0.97 [20].

Zoom In Zoom Out Reset image size

At the high-mass end, the overlap is dominated by NR data, and as in figure 7, the overlaps are high without needing to move off the signal mass. At the low-mass end, the same result would be expected in a pure pN/pN comparison although there is enough of the hybridization in-band to reduce the overlaps. However, changing the mass introduces a phase difference that accumulates over all the cycles in-band, and so higher overlaps cannot be achieved. The result is optimal mass value close to the correct mass value, but with a lower overlap. In the middle region, these factors compete. At higher masses, the overlap is reduced less by changing the mass and so the recovered value of the mass can stray further from the injected value. As the hybridization passes out of band, this adjustment is no longer needed. The same general behavior can be seen in comparisons between non-spinning, unequal-mass (q = 2) waveforms and equal-mass spinning waveforms, as shown in figure 11, which shows overlaps between waveforms with identical parameters and hybridized with the same pN approximants maximized over the mass of one waveform. We can infer from these results that although we have only made a crude estimate of the parameter bias, the accuracy of the waveforms (excluding the uncertainty in the pN approximants) is extremely high. There are many factors that may bias parameter estimation of real signals, including uncertainties in the detector calibration, noise in the detectors and errors in pN waveforms used as templates. Subsequent NINJA-2 data-analysis studies will attempt to quantify the effects of these factors.

Zoom In Zoom Out Reset image size

The efficiency of searches for GW signals from BH binaries and the accuracy with which the source parameters are estimated will depend crucially on progress in incorporating information from approximation methods and NR into the waveform models used in data-analysis algorithms. Even in the non-spinning case, a recent extensive comparison of different pN approximants [133] concludes that for masses above 12M_⊙ NR simulations of the last orbits and merger are required for the construction of optimal detection templates. This need is expected to be even greater when spinning binaries are included. And of course, NR simulations will be yet more important for parameter estimation.

Injections of inspiral–merger–ringdown (IMR) waveforms from analytical waveform models are already used to calibrate the analysis of detector data [136]. However, the synthesis of NR results into waveform models typically lags several years behind the most complete current sets of NR waveforms. The NINJA-2 project will consider direct injections of hybrid pN–NR waveforms, which will allow the use of waveforms that have not (yet) been used in constructing analytical models and will avoid any additional modeling errors. Using the best available IMR waveforms will be particularly important for parameter estimation.

The first goal of the second NINJA project has been to review the submitted hybrid waveforms. In this paper, we have summarized the requirements that the submitted waveforms must meet (see section 1). All submitted waveforms meet these requirements, within some minor caveats that are detailed in section 1. We have also demonstrated the validity of the hybrid waveforms for use in the NINJA-2 project with more detailed consistency checks of the submitted data.

The validation of and comparisons between the submissions, which we have presented here, are a major feature of NINJA-2 that was not present in NINJA-1 and will be indispensable when analyzing the results from systematic injection studies over the course of the NINJA-2 project. A total of 56 waveforms from eight different groups were contributed to NINJA-2, corresponding to 40 distinct NR waveforms, and 29 different configurations of mass ratio and spins. For six configurations, multiple numerical waveforms were submitted, and 16 numerical waveforms were hybridized with two different pN approximants (TaylorT1 and TaylorT4). This has allowed waveform comparisons and tests of the accuracy standards, as discussed in detail in section 5, and summarized below. For each submission, significant further information has been included in metadata files as specified in [61], e.g., hybridization frequencies, eccentricities, or literature references. This allowed us to automatically generate information, such as table 1, and has been proved crucial in systematically analyzing the data set. Verifying these submissions, preparing a consistent data set and evaluating their accuracy by comparing submissions with the same mass ratio and spin, has been a formidable task, and in this paper, we only report on the (ℓ, m) = (2, 2) spherical-harmonic modes. Based on our various checks and comparisons as reported in this paper, we judge the submitted waveforms suitable for the NINJA-2 project. Parameter estimation places particular stringent demands on waveform templates. We believe the NINJA-2 waveforms to be of sufficient quality to estimate size and shape of maximum likelihood contours. However, the large truncation error of the pN expansion in the cycles before and during hybridization will induce systematic errors; further investigations will be necessary before these waveforms can be used to check the accuracy of the estimated parameters themselves (i.e. the location of the maximum likelihood contour, rather than its shape).

As shown in figure 2, the submissions primarily cover only two lines in the (q, χ) plane, leaving large regions of parameter space unexplored. We also do not consider any precessing signals, which adds several dimensions to the parameter space. These issues will be explored in future NINJA projects. Extending our analysis to non-dominant spherical-harmonic modes, and to a larger volume of parameter space, will mark a significant challenge for the future, and will require a significant extension of our methods of analysis, and of automatizing them, and will constitute a substantial research project by itself.

An important shortcoming in our analysis of waveform accuracy is that it is currently not possible to accurately estimate the pN truncation error. A reasonable approach seems to be the comparison of pN approximants, which are accurate to the highest known order, but differ in terms beyond this order (3.5 pN without spins). These deviations still depend strongly on the choice of pN approximants used in any comparison, and on the mass ratio and spins. For example, in figures 6, 7 and 10, we compare the use of the TaylorT1 and TaylorT4 approximants for the mass ratio q = 1 and zero spins; figure 7 also includes q ≠ 1. Previous work implies that TaylorT4 performs exceptionally well in these cases [25, 27], apparently by coincidence, but beyond that it is not possible to make precise quantitative statements; this figure simply illustrates the general level of mismatch error that we may expect between various pN approximants across a range of binary masses.

We caution that the comparisons presented in section 5 were made only for subset of the black-hole configurations, which do not include the most extreme cases and not the cases with unequal Kerr parameters. Further work will be necessary to establish with similar confidence that all NINJA-2 hybrid waveforms are of similar quality.

Our waveform comparisons in section 5 are consistent with previous results: mismatches between waveforms are dominated by the choice of pN approximant, while NR and hybridization errors are far smaller. Hybridization choices and methods certainly affect overlaps, although the mismatch due to hybridization appears to be less than a fraction of a percent. This is not expected to have any noticeable effect on detection, but is likely to impact parameter estimation. The degree to which these effects will bias searches is one of the key question we hope NINJA-2 will be able to answer. In particular, theoretical studies such as the ones presented in this paper cannot replace a complete parameter-estimation analysis, and it is not clear how theoretical studies based on Gaussian noise can predict the performance of GW searches in real noise. Another important goal of the NINJA-2 project is thus to build up experience in comparing simple theoretical studies with the full GW-search plus parameter-estimation pipeline exercised in actual observations of GW events.

We thank Ilya Mandel for helpful discussions. We gratefully acknowledge support from the National Science Foundation under NSF grants PHY-1040231, PHY-0600953, PHY-0847611, AST-1028087, DRL-1136221, OCI-0832606, PHY-0903782, PHY-0929114, PHY-0969855, AST-1002667, PHY-0650377, PHY-0963136, PHY-0855315, PHY-0969111, PHY-1005426, PHY-0601459, PHY-1068881, PHY-1005655, PHY-0653550, PHY-0955773, PHY-0653443, PHY-0855892, PHY-0914553, PHY-0941417, PHY-0903973, PHY-0955825, NSF cooperative agreement PHY-0757058, by NASA grants 07-ATFP07-0158, NNX07AG96G, NNX10AI73G, NNX09AF96G, NNX09AF97G, by Marie Curie Grants of the 7th European Community Framework Programme FP7-PEOPLE-2011-CIG CBHEO number 293412, by the DyBHo-256667 ERC Starting grant, and MIRG-CT-2007-205005/PHY, and Science and Technology Facilities Council grants ST/H008438/1 and ST/I001085/1. Further funding was provided by the Sherman Fairchild Foundation, NSERC of Canada, the Canada Research Chairs Program, the Canadian Institute for Advanced Research, Govern de les Illes Balears, the Ramón y Cajal Programme of the Ministry of Education and Science of Spain, contracts AYA2010-15709, CSD2007-00042, FIS2011-30145-C03-03, CSD2009-00064 and FPA2010-16495 of the Spanish Ministry of Science and Innovation, grant AGAUR-2009-SGR-935, the Royal Society, the German Research Foundation, grant SFB/Transregio 7, the German Aerospace Center for LISA Germany, and the Research Corporation for Science Advancement. Computations were carried out on Teragrid machines Lonestar, Ranger, Trestles and Kraken under Teragrid allocations TG-PHY060027N, TG-MCA99S008, TG-PHY090095, TG-PHY100051, TG-PHY990007N, TG-PHY090003 and TG-MCA08X009. Computations were also performed on the clusters "HLRB-2' at LRZ Munich, 'NewHorizons' at RIT (funded by NSF grant numbers AST-1028087, DMS-0820923 and PHY-0722703), 'Zwicky' at Caltech (funded by NSF MRI award PHY-0960291), 'Finis Terrae' (funded by CESGA-ICTS-2010-200 and 221), 'Caesaraugusta' (funded by BSC grant numbers AECT-2011-2-0006, AECT-2011-3-0007, AECT-2012-1-0008), 'MareNostrum' (funded by BSC grant numbers AECT-2009-2-0017, AECT-2010-1-0008, AECT-2010-2-0013, AECT-2010-3-0010, AECT-2011-1-0015, AECT-2011-2-0012), 'VSC' in Vienna (funded by FWF grant P22498), 'Force' at GaTech, and on the GPC supercomputer at the SciNet HPC Consortium [137]; SciNet is funded by the Canada Foundation for Innovation under the auspices of Compute Canada; the Government of Ontario; Ontario Research Fund, Research Excellence; and the University of Toronto.

The accuracy of hybrid waveforms depends very sensitively on the accuracy of the pN waveforms with which they are constructed. For this reason, the NINJA-2 collaboration invested significant effort into ensuring that all contributions used the most current pN information available. Mathematica notebooks containing the full expressions and derivations of the various approximants are provided in the ancillary files available with this paper online. Here, we describe the techniques.

The pN approximants used in this paper solve for the orbital motion to high accuracy by assuming that the motion is an adiabatic quasicircular inspiral and by assuming that the loss of orbital binding energy E during inspiral is balanced by the flux of energy in gravitational radiation to infinity $\mathcal {F}$ and the flux of energy going into the individual black holes caused by tidal heating $\dot{M}$. The first assumption means that the motion can be described completely by the orbital phase function Φ(t). We then define the pN expansion parameter v := (M dΦ/dt)^1/3, where M is the sum of the apparent-horizon masses of the black holes. The orbital binding energy, gravitational-wave flux and tidal heating can be expressed as functions of this parameter. Thus, the energy-balance equation can be written as $\dot{E} + \mathcal {F} + \dot{M} = 0$. Using the chain rule to rewrite $\dot{E}$, we can rearrange this as an expression for dv/dt and include the expression for dΦ/dt to obtain a complete system of ordinary differential equations describing the motion of the binary:

$$\begin{equation} \frac{\mathrm{d}v}{\mathrm{d}t} = - \frac{\mathcal {F}(v) + \dot{M}(v)}{E^{\prime }(v)}, \qquad \frac{\mathrm{d}\Phi }{\mathrm{d}t} = \frac{v^{3}}{M}. \end{equation} \tag{ A.1 }$$

The formulation of the pN approximants, then, comes down to writing down the expressions for E(v), $\mathcal {F}(v)$ and $\dot{M}(v)$ and integrating the system for v(t) and Φ(t).

At leading order, the expressions for E and $\mathcal {F}$ are

$$\begin{equation} E(v) = -\frac{M \eta \, v^2}{2} \qquad {\rm and} \qquad \mathcal {F}(v) = \frac{32}{5}\, v^{10}\, \eta ^2, \end{equation} \tag{ A.2 }$$

where higher order terms include additional factors of v. The additional terms are currently known up to v⁷ (3.5 pN) for non-spinning systems and to lower order for spinning systems. The tidal heating $\dot{M}$ is equivalent to a 2.5 pN term in the flux expression. The full expressions for E, $\mathcal {F}$ and $\dot{M}$ are calculated in [138–144], with a collection of errata in [61], and are given explicitly in the accompanying Mathematica notebook PNOrbitalPhase.nb, which can be found in the ancillary material with this paper online. For all the results in this paper, the full expressions were used.

To integrate the system of ordinary differential equations, various methods have been developed, each of which should be equivalent at the level of accuracy of the pN approximation. These 'approximants' have been given the names TaylorT1 to TaylorT4 [25, 145, 146]. The TaylorT1 approximant is constructed by directly evaluating the expressions in equation (A.1), which are then integrated numerically. For the TaylorT4 approximant, the right-hand side of the expression for dv/dt is first re-expanded in a Taylor series and truncated at the appropriate 3.5 pN order, which is then integrated numerically. The difference between these approximants is at the level of the 4 pN uncertainty.

Alternatively, we can find analytical formulas, rather than integrating numerically. The TaylorT2 approximant is derived by taking the (multiplicative) inverse of the first expression in equation (A.1) to construct a new system:

$$\begin{equation} \frac{\mathrm{d}t}{\mathrm{d}v} = - \frac{E^{\prime }(v)}{\mathcal {F}(v) + \dot{M}(v)}, \qquad \frac{\mathrm{d}\Phi }{\mathrm{d}v} = \frac{\mathrm{d}\Phi }{\mathrm{d}t} \frac{\mathrm{d}t}{\mathrm{d}v} = - \frac{v^{3}}{M}\, \frac{E^{\prime }(v)}{\mathcal {F}(v) + \dot{M}(v)}. \end{equation} \tag{ A.3 }$$

We re-expand the right-hand sides of these expressions in the Taylor series, truncate at the appropriate 3.5 pN order and integrate with respect to v, obtaining expressions for t(v) and Φ(v). This parametrically determines the phase of the binary as a function of time, with v being the independent parameter. Finally, the TaylorT3 approximant is constructed by inverting the series t(v) to obtain v(t). There is a 3 pN logarithmic term in t(v), which must be treated as a constant in order to invert the series. Once this is done, the series for v(t) can be inserted into dΦ/dt = v³/M, which can be integrated analytically to find Φ(t).

Using any of these approximants, the resulting orbital phase and frequency allow us to calculate the metric perturbation function h. The perturbation falls off at leading order as 1/r and varies with angle, typically being decomposed in spin-weight s = −2 spherical harmonics [61]. The dominant mode of this decomposition is generally the (ℓ, m) = (2, 2) mode. At leading order, we have²⁰

$$\begin{equation} h_{2,2}(v, \Phi ) = \frac{M}{r}\, \sqrt{\frac{\pi }{5}}\, 8\, \eta \, v^{2}\, \mathrm{e}^{-\mathrm{i}2 \Phi }, \end{equation} \tag{ A.4 }$$

where higher order terms include additional factors of v. The additional terms are currently known up to v⁶ (3 pN) for non-spinning systems [147] and to lower order for spinning systems [61, 148]. The full expressions used in this paper, including other modes, are given in the accompanying Mathematica notebook PNWaveform.nb found in the ancillary materials with this paper online.

The approximants just discussed describe the pN waveform in the time domain. For many purposes, it can be useful to have expressions in the frequency domain. The frequency-domain approximant TaylorF2 is obtained from TaylorT2 using the stationary-phase approximation (SPA) [149], together with the approximation that the gravitational-wave frequency is just twice the orbital frequency, so f = v³/π M. Then, the frequency-domain waveform is given by

$$\begin{equation} \fl \tilde{h}_{\ell ,m}(f) = h_{\ell ,m}(v, 0)\, \sqrt{ \frac{2\pi M}{3 m v^{2}\, \dot{v}} }\, \exp \lbrace \mathrm{i}[ 2 v^{3}\, t(v)/M - m\Phi (v) -\pi /4 ] \rbrace , \end{equation} \tag{ A.5 }$$

where $\dot{v}$ is given in equation (A.1), and t(v) and Φ(v) are the results of the TaylorT2 approximant.