Targeted numerical simulations of binary black holes for GW170104
Abstract
In response to LIGO’s observation of GW170104, we performed a series of full numerical simulations of binary black holes, each designed to replicate likely realizations of its dynamics and radiation. These simulations have been performed at multiple resolutions and with two independent techniques to solve Einstein’s equations. For the nonprecessing and precessing simulations, we demonstrate the two techniques agree mode by mode, at a precision substantially in excess of statistical uncertainties in current LIGO’s observations. Conversely, we demonstrate our full numerical solutions contain information which is not accurately captured with the approximate phenomenological models commonly used to infer compact binary parameters. To quantify the impact of these differences on parameter inference for GW170104 specifically, we compare the predictions of our simulations and these approximate models to LIGO’s observations of GW170104.
I Introduction
The LIGOVirgo Collaboration (LVC) has reported the confident discovery of three binary black hole (BBH) mergers via gravitational wave (GW) radiation: GW150914Abbott et al. (2016a) and GW151226Abbott et al. (2016b) from the first observing run O1Abbott et al. (2016c), and GW170104 Abbott et al. (2017), GW170608 Abbott et al. (2017a), and GW170814 Abbott et al. (2017b) from the second observing run. The parameters of these detections were inferred by comparing the data to stateoftheart approximate models Taracchini et al. (2014); Pürrer (2014); Hannam et al. (2014). A reanalysis of GW150914 implementing full numerical relativity (NR)Abbott et al. (The LIGO Scientific Collaboration and the Virgo Collaboration) (2016) simulations helped to better constrain the mass ratio of the system. This is due to the fact that NR waveforms include physics omitted by current approximate models, notably higher order modes and accurate precession effects. A full description of this methodology, including detailed tests of systematic errors and parameter estimation improvements, can be found in Lange et al. (2017).
This paper is organized as follows. In Section II, we describe the two independent techniques we use to solve Einstein’s equations numerically for the evolution of binary black hole spacetimes. In Section III, we describe the binary’s parameters selected for detailed followup, our simulations of these proposed initial conditions, and detailed comparisons between our paired results, for both nonprecessing and precessing simulations. We also contrast our simulations’ radiation with the corresponding results derived from the approximate phenomenological models used by LIGO for parameter inference. In Section IV, we directly compare our simulations to GW170104. These comparisons provide both a scalar measure of how well each simulation agrees with the data (a marginalized likelihood), as well as the bestfitting reconstructed waveform in each instrument Abbott et al. (The LIGO Scientific Collaboration and the Virgo Collaboration) (2016); Lange et al. (2017). Using our reconstructed waveforms, we graphically demonstrate that our simulations agree with each other and the data, with simulation differences far smaller than the residual noise in each instrument. Using these real observations as a benchmark for model quality, we then quantify how effectively our simulations reproduce the data, compared to the results of approximate and phenomenological models at the same parameters. Since our simulations parameters were selected using these approximate and phenomenological models, we also have the opportunity to assess how effectively they identified the optimal binary parameters. In Section V we discuss the prospects for future targeted simulations in followup of LIGO observations.
Ii Full Numerical Evolutions
The breakthroughs Pretorius (2005); Campanelli et al. (2006); Baker et al. (2006) in numerical relativity allowed for detailed predictions for the gravitational waves from the late inspiral, plunge, merger and ringdown of black hole binary systems. Catalogs of the simulated waveforms are publicly available Mroue et al. (2013); Jani et al. (2016); Healy et al. (2017a) for its use for BBH parameter estimation Lange et al. (2017), as well as for determining how the individual masses and spins of the orbiting binary relate to the properties of the final remnant black hole produced after merger. This relationship Healy et al. (2014) can be used as a consistency check for the observations of the inspiral and, independently, the mergerringdown signals as tests of general relativity Ghosh et al. (2016); Abbott et al. (2016d, c).
ii.1 Simulations using finitedifference, movingpuncture methods
In order to make systematic studies and build a data bank of full numerical simulations, e.g., Healy et al. (2017a), it is crucial to develop efficient numerical algorithms, since large computational resources are required. The Rochester Institute of Technology (RIT) group evolved the BBH data sets described below using the LazEv Zlochower et al. (2005) implementation of the moving puncture approach Campanelli et al. (2006); Baker et al. (2006) with the conformal function suggested by Ref. Marronetti et al. (2008). For those runs, they used centered, sixthorder finite differencing in space Lousto and Zlochower (2008) and a fourthorder Runge Kutta time integrator (the code does not upwind the advection terms) and a 5thorder KreissOliger dissipation operator.
The LazEv code uses the EinsteinToolkit Löffler et al. (2012); ein / Cactus cac / Carpet Schnetter et al. (2004) infrastructure. The Carpet mesh refinement driver provides a “moving boxes” style of mesh refinement. In this approach, refined grids of fixed size are arranged about the coordinate centers of both holes. The Carpet code then moves these fine grids about the computational domain by following the trajectories of the two BHs.
To compute the initial low eccentricity orbital parameters RIT used the postNewtonian techniques described in Healy et al. (2017b) and then generated the initial data based on these parameters using approach Brandt and Brügmann (1997) along with the TwoPunctures Ansorg et al. (2004) code implementation.
The LazEv code uses AHFinderDirect Thornburg (2004) to locate apparent horizons, and measures the magnitude of the horizon spin using the isolated horizon (IH) algorithm detailed in Ref. Dreyer et al. (2003) and as implemented in Ref. Campanelli et al. (2007). The horizon mass is calculated via the Christodoulou formula where , is the surface area of the horizon, and is the spin angular momentum of the BH (in units of ).
The radiated energy, linear momentum, and angular momentum, were measured in terms of the radiative Weyl Scalar , using the formulas provided in Refs. Campanelli and Lousto (1999); Lousto and Zlochower (2007), Eqs. (22)(24) and (27) respectively. However, rather than using the full , it was decomposed it into and modes and dropping terms with . The formulas in Refs. Campanelli and Lousto (1999); Lousto and Zlochower (2007) are valid at . To obtain the waveform and radiation quantities at infinity, the perturbative extrapolation described in Ref. Nakano et al. (2015) was used.
For the RIT simulations, different resolutions are denoted by NXXX where XXX is either 100, 118, or 140 for low, medium, and high resolutions, respectively. This number is directly related to the wavezone resolution in the simulation. For instance, N100 has a resolution of in the wavezone (where observer extraction takes place, preliminary to perturbative extrapolation to infinity via Nakano et al. (2015)), and N140 has . In each case UID#15, there are 10 levels of refinement in all and the grids followed a pattern close to those described in Healy et al. (2017c).
Other groups using the moving punctures Campanelli et al. (2006); Baker et al. (2006) formalism with finite difference methods are Georgia Institute of Technology (GT) Jani et al. (2016) and those based on BAM Brügmann et al. (2008). The GT Jani et al. (2016) simulations were obtained with the Maya code Herrmann et al. (2007a, b); Hinder et al. (2008); Healy et al. (2009a); Hinder et al. (2010); Healy et al. (2009b, 2010); Bode et al. (2010), which is also based on the BSSN formulation with moving punctures. The grid structure for each run consisted of 10 levels of refinement provided by Carpet Schnetter et al. (2004), a mesh refinement package for Cactus cac . Each successive level’s resolution decreased by a factor of 2. Sixthorder spatial finite differencing was used with the BSSN equations implemented with Kranc Husa et al. (2006).
ii.2 Simulations using pseudospectral, excision methods
Simulations labeled SXS are carried out using the Spectral Einstein Code (SpEC) SpE used by the Simulating eXtreme Spacetimes Collaboration (SXS). Given initial BBH parameters, a corresponding weighted superposition of two boosted, spinning KerrSchild black holes Lovelace et al. (2008) is constructed, and then the constraints are solved York (1999); Pfeiffer et al. (2003); Ossokine et al. (2015) by a pseudospectral method to yield quasiequilibrium Caudill et al. (2006); Lovelace et al. (2008) initial data. Small adjustments in the initial orbital trajectory are made iteratively to produce initial data with low eccentricity Pfeiffer et al. (2007); Buonanno et al. (2011); Mroué and Pfeiffer (2012).
The initial data are evolved using a firstorder representation Lindblom et al. (2006) of a generalized harmonic formulation Friedrich (1985); Garfinkle (2002); Pretorius (2005) of Einstein’s equations, and using damped harmonic gauge Lindblom and Szilágyi (2009); Choptuik and Pretorius (2010); Szilágyi et al. (2009). The equations are solved pseudospectrally on an adaptivelyrefined Lovelace et al. (2011); Szilágyi (2014) spatial grid that extends from pureoutflow excision boundaries just inside apparent horizons M. A. Scheel, M. Boyle, T. Chu, L. E. Kidder, K. D. Matthews and H. P. Pfeiffer (2009); Szilágyi et al. (2009); Hemberger et al. (2013); Ossokine et al. (2013); Scheel et al. (2015) to an artificial outer boundary. Adaptive timestepping automatically achieves time steps of approximately the Courant limit.
On the Cal State Fullerton cluster, ORCA, the simulation achieved a typical evolution speed of for the highest resolution (here we measure simulation time in units of , the total mass of the binary). After the holes merge, all variables are automatically interpolated onto a new grid with a single excision boundary inside the common apparent horizon M. A. Scheel, M. Boyle, T. Chu, L. E. Kidder, K. D. Matthews and H. P. Pfeiffer (2009); Hemberger et al. (2013), and the evolution is continued. Constraintpreserving boundary conditions Lindblom et al. (2006); Rinne (2006); Rinne et al. (2007) are imposed on the outer boundary, and no boundary conditions are required or imposed on the excision boundaries.
We use a pseudospectral fastflow algorithm Gundlach (1998) to find apparent horizons, and we compute spins on these apparent horizons using the approximate Killing vector formalism of Cook, Whiting, and Owen Cook and Whiting (2007); Owen (2007).
Gravitational wave extraction is done by three independent methods: direct extraction of the NewmanPenrose quantity at finite radius M. A. Scheel, M. Boyle, T. Chu, L. E. Kidder, K. D. Matthews and H. P. Pfeiffer (2009); Boyle et al. (2007); Pfeiffer et al. (2007), extraction of the strain by matching to solutions of the ReggeWheelerZerilliMoncrief equations at finite radius Buchman and Sarbach (2007); Rinne et al. (2009), and CauchyCharacteristic Extraction Bishop et al. (1997); Winicour (2005); Gomez et al. (2007); Reisswig et al. (2013); Handmer and Szilágyi (2015). The latter method directly provides gravitational waveforms at future null infinity, while for the former two methods the waveforms are computed at a series of finite radii and then extrapolated to infinity Boyle and Mroué (2009). Differences between the different methods, and differences in extrapolation algorithms, can be used to as error estimates on waveform extraction. These waveform extraction errors are important for the overall error budget of the simulations, and are typically on the order of, or slightly larger than, the numerical truncation error Taylor et al. (2013); Chu et al. (2016). In this paper, the waveforms we compare use ReggeWheelerZerilliMoncrief extraction and extrapolation to infinity. We have verified that our choice of extrapolation order does not significantly affect our results. We have also checked that corrections to the wave modes Boyle (2016a) to account for a small drift in the coordinates of the center of mass have a negligible effect on our results.
Iii Simulations of GW170104
We extracted the maximum a posteriori (MaP) parameters from (preliminary) Bayesian posterior inferences performed by the LIGO Scientific Collaboration and the Virgo Collaboration, using different waveform models Abbott et al. (2017c); Veitch et al. (2015). As described in Appendix B, this point parameter estimate is one of a few wellmotivated and somewhat different choices for followup parameters; however, as described in that appendix, we estimate that the specific choice we adopt will not significantly change our principal results. Table 1 shows parameters simulated with numerical relativity.
Run  [Hz]  Approximant  

#1  58.49  24  0.8514  ( 0, 0, 0.7343)  ( 0, 0, 0.8278)  SEOBNRv4ROM 
#2  58.72  24  0.5246  ( 0.1607, 0.1023, 0.0529 )  ( 0.3623, 0.5679, 0.3474 )  SEOBNRv3 
#3 (*)  62.13  20  0.4850  ( 0.0835, 0.4013, 0.3036 )  ( 0.3813, 0.7479, 0.1021 )  IMRPhenomPv2 
#4  53.46  20  0.7147  ( 0, 0, 0.2205)  ( 0, 0, 0.7110)  SEOBNRv4ROM 
#5  59.11  20  0.4300  ( 0, 0, 0.3634)  ( 0, 0, 0.1256)  IMRPhenomD 
Spin Conventions: are specified in a frame where (i) , i.e. the Newtonian orbital angular momentum is along the zaxis. (ii) the vector pointing from to is the xaxis, (1,0,0). Note that the orientation of is essentially undetermined by parameter estimation (PE) methods, so the choice (ii) is meant to break this degeneracy to arrive at concrete parameters. In other words, the spincomponents given below are those consistent with Eqs. (9)(11) of
^{1}^{1}1https://dcc.ligo.org/DocDB/0123/T1500606/006/
NRInjectionInfrastructure.pdf
The label, UID#15, of the simulation identifies which parameters we are using in following up as given by the initial data from Table 1. For aligned spin runs, where spinvectors are preserved, the initial orbital frequency may be smaller than . For precessing runs, because we target a certain spin configuration at , this forces the initial NRfrequency to be identical (or only somewhat smaller than) the reference frequency. is the initial orbital separation of the NR run in geometric units and the mass ratio and intrinsic spins are denoted by . Due to the way NR simulations are setup the initial parameters can change due to the presence of junk radiation and/or imperfections in setting up initial data, therefore these quantities should be reported as afterjunk masses/spins, ideally extracted at the reference frequency. For precessing runs, in particular, the spincomponents should be specified at the reference frequency, following the convention , with and computed at the reference frequency, too. We also provide , the orbital eccentricity. For instance, the actual initial data as measured for the RIT’s followup simulations are described in Table 2
Run  

#1  7.9168  6.7407  2.829e4  0.07467  0.3196  0.3056  0.4599  0.5401  0.9928  0.7343  0.8267  6e4 
#2  7.8211  4.1029  4.837e4  0.07837  0.3277  0.4400  0.3441  0.6559  0.9922  0.1977  0.7580  1e3 
#3  8.7720  4.2543  3.316e4  0.07160  0.2796  0.3584  0.3266  0.6734  0.9930  0.5101  0.8445  1e3 
#4  8.4742  6.0567  2.918e4  0.07421  0.3991  0.4219  0.4168  0.5832  0.9928  0.2205  0.7110  2e4 
#5  9.4395  4.0593  2.718e4  0.06698  0.2753  0.6865  0.3007  0.6993  0.9933  0.3634  0.1256  5e4 
UID1N100  0.459757  0.539780  0.851718  (0,0,0.737493)  (0,0,0.829064)  0.002586 

UID1N118  0.459758  0.539801  0.851718  (0,0,0.737473)  (0,0,0.829040)  0.002587 
UID1N140  0.459758  0.539801  0.851718  (0,0,0.737464)  (0,0,0.829030)  0.002587 
UID1L2  0.459913  0.540183  0.851401  (0,0,0.734254)  (0,0,0.827442)  0.003135 
UID1L3  0.459902  0.540171  0.851400  (0,0,0.734138)  (0,0,0.827657)  0.003139 
UID1L4  0.459907  0.540176  0.851402  (0,0,0.734106)  (0,0,0.827690)  0.003139 
UID2N100  0.344090  0.655693  0.524773  ( 0.119494, 0.144676, 0.063539 )  ( 0.675724, 0.181910, 0.293427 )  0.004010 
UID2N118  0.344091  0.655693  0.524774  ( 0.118992, 0.144981, 0.063537 )  ( 0.676089, 0.180798, 0.293141 )  0.004010 
UID2N140  0.344091  0.655694  0.524774  ( 0.118754, 0.145124, 0.063541 )  ( 0.676262, 0.180273, 0.293005 )  0.004010 
UID2L1  0.344082  0.655876  0.524614  ( 0.101171, 0.156396, 0.066135 )  ( 0.693687, 0.096203, 0.290389 )  0.004001 
UID2L2  0.344069  0.655966  0.524523  ( 0.100390, 0.156863, 0.066508 )  ( 0.693803, 0.095200, 0.289578 )  0.004013 
UID2L3  0.344072  0.655958  0.524534  ( 0.100596, 0.156736, 0.066591 )  ( 0.693727, 0.095442, 0.289439 )  0.004010 
UID3N100  0.326605  0.672963  0.485325  ( 0.360515, 0.132028, 0.338700 )  ( 0.685165, 0.491281, 0.068777 )  0.003308 
UID3N118  0.326606  0.672963  0.485326  ( 0.360458, 0.131982, 0.338637 )  ( 0.685154, 0.491237, 0.068765 )  0.003308 
UID3N140  0.326607  0.672964  0.485326  ( 0.360429, 0.132025, 0.338591 )  ( 0.685106, 0.491276, 0.068763 )  0.003308 
UID3L1  0.326576  0.673382  0.484978  ( 0.358221, 0.095241, 0.350438 )  ( 0.749729, 0.387741, 0.052141 )  0.003330 
UID3L2  0.326598  0.673451  0.484962  ( 0.358434, 0.094398, 0.350478 )  ( 0.750204, 0.386567, 0.051977 )  0.003339 
UID3L3  0.326583  0.673462  0.484931  ( 0.352293, 0.114550, 0.350800 )  ( 0.727233, 0.427282, 0.059549 )  0.003308 
UID4N100  0.416817  0.583036  0.714908  (0,0,0.221608)  (0,0,0.712189)  0.002629 
UID4N118  0.416819  0.583037  0.714911  (0,0,0.221602)  (0,0,0.712184)  0.002630 
UID4N140  0.416820  0.583037  0.714912  (0,0,0.221600)  (0,0,0.712182)  0.002630 
UID4L1  0.416790  0.583195  0.714666  (0,0,0.220660)  (0,0,0.710930)  0.002635 
UID4L2  0.416809  0.583214  0.714675  (0,0,0.220453)  (0,0,0.710967)  0.002635 
UID4L3  0.416817  0.583203  0.714703  (0,0,0.220427)  (0,0,0.710930)  0.002619 
UID5N100  0.300721  0.699282  0.430043  (0,0,0.366067)  (0,0,0.125360)  0.002932 
UID5N118  0.300722  0.699284  0.430043  (0,0,0.366031)  (0,0,0.125352)  0.002932 
UID5N140  0.300723  0.699285  0.430043  (0,0,0.366015)  (0,0,0.125348)  0.002932 
UID5L1  0.300697  0.699261  0.430021  (0,0,0.363371)  (0,0,0.125641)  0.002900 
UID5L2  0.300708  0.699223  0.430061  (0,0,0.363424)  (0,0,0.125550)  0.002864 
UID5L3  0.300722  0.699178  0.430107  (0,0,0.363345)  (0,0,0.125616)  0.002855 
For followup #1, the initial spurious burst contains a nonnegligible kick which causes a center of mass drift of approximately 0.65M over the 4000M of evolution. Because of this, information from the dominant modes leak into the other modes, particularly the odd modes. To reduce this effect, we can recalculate the modes by finding the average rest frame of the binary. We calculate the average velocity of the center of mass of the binary (from ) over the inspiral and then boost the waveform in the opposite direction. This is done using Eqs. (45) in Boyle (2016b) and Eqs. (78) in Kelly and Baker (2013). Note that this does not change the physical waveforms, only how they are spread over modes.
For the RIT simulations, the initial data parameters in Table 2 for the nonprecessing systems 1, 4, and 5 were determined by choosing the starting frequency just below the reference frequency. This gives the gauge time to settle, and since the spins do not change, this gives us a cleaner waveform once we hit the reference frequency. For the precessing simulations 2 and 3, since the spins will now evolve, we determine the initial data parameters by choosing the initial spins at the specified reference frequency. The initial data used by SXS, being not conformally flat have less spurious radiation content than the BowenYork data and hence produce a different set of masses and spins after settling down. See Table 3 for the specific values of each simulation by the two kind of initial data families. This process can be iterated to get closer initial parameters for each approach, although it requires some extra evolution time and coordination to reach a fractional agreement below . This process has been followed for UID#1, but not for the other cases, in particular the two precessing ones #2 and #3, and hence the differences, for instance, displayed in Fig. 1 for the precessing case #3.
The SXS simulations used in this work have been assigned SXS catalog numbers BBH:0626 (UID1), BBH:0627 (UID2), BBH:0628 (UID3), BBH:0625 (UID4), and BBH:0631 (UID5).
Each simulation has an asymptotic frame relative to which we extract . In all cases used here, this axis corresponds to the () axis of the simulation frame. For all simulations, this axis also agrees with the orbital angular momentum axis at the start of the evolution.
iii.1 Outgoing radiation very similar for different NR methods
Following previous studies Lovelace et al. (2016), we compare the outgoing radiation mode by mode, using an observationallydriven measure: the overlap or match. The black and grey lines in Figures 2 and 3 show the match between the two simulations’ (RITSXS and RITGT respectively) (2,2) modes, as a function of the minimum frequency used in the match. In this calculation, we use a detector noise power spectrum appropriate to GW170104, and a total mass as given in Table 1. By increasing the minimum frequency, we increasingly omit the earliest times in the signal, first eliminating transient startup effects associated due to finite duration and eventually comparing principally the merger signals from the two black holes. For comparison, the red, blue, and yellow lines show the corresponding matches between RIT, SXS, and GT simulations respectively and effective one body models with identical parameters (faithfulness study). In Figure 2, which illustrates only nonprecessing simulations, these comparisons are made to the nonprecessing model SEOBNRv4 Bohé et al. (2017). In Figure 3, which targets the two precessing UIDs, we instead compare to SEOBNRv3, which approximates some precession effects. For both nonprecessing and precessing simulations, these figures show that the different NR groups’ simulations produce similar radiation, with mismatches even at the longest durations considered. By contrast, comparisons with SEOBNRv4 and SEOBNRv3 show that these models do not replicate our simulations’ results, particularly for precessing binaries.
To demonstrate good agreement beyond the (2,2) mode for precessing simulations, for multiple resolutions, Tables 4 and 5 systematically compare all modes between RIT and SXS. The match calculations in this Table are performed using a strain noise power spectral densities (PSD) characterizing data near GW170104. Following Lovelace et al. (2016), one phase and timeshift is computed by maximizing the overlap of the (2,2) mode; this phase and timeshift is then applied to all other modes without any further maximization. Table 4 shows a resolution test: the match between RIT and SXS simulations, as a function of RIT simulation resolution. As the most challenging precessing case, UID3 is shown by default Except for the modes, all the simulations show good agreement modebymode, for all resolutions.
Based on Figures 2, 3 and Table 4 we anticipate the lowest productionquality NR resolution (N100 for RIT; L3 for SXS) will usually be sufficient to go well beyond the accuracy of approximate and phenomenological models. To elaborate on this hypothesis, Table 5 shows the modebymode overlaps between these two lowest NR resolutions. The modes agree without exception. Good agreement also exists for the most significant modes up to . On the other hand, the last columns of Table 4 shows that the rough agreement between NR and models for the modes (2,2) displayed in Fig. 3, notably worsens when looking at other than the leading modes.
In Table 6, we also provide a comparison of the remnant properties, ie. final mass, spin and recoil velocity of the final, merged, black hole, as computed by the two NR methods and for a set of three increasing resolutions. We observe good agreement and convergence of their values. In the case of the RIT runs, a nearly 4th order convergence with resolution for the recoil velocity of the remnant is displayed. The same convergence properties are shown for the final mass and spin, despite these quantities being overresolved at these resolutions. Those remnant quantities are also important to model fitting formulae Healy and Lousto (2017); JimńezForteza et al. (2017) to be used to infer the final black hole properties from the binary parameters and thus serve as a test of the general theory of relativity, as in Ghosh et al. (2016). The phenomenological approximate waveform models Bohé et al. (2017); Blackman et al. (2017a); London et al. (2017) also benefit from information of the final remnant properties as one of their inputs and can thus produce a more accurate precessing and including higher modes model.
In conclusion we see that the typical production NR simulations are well into the convergence regime and produce accurate enough waveforms for all practical applications of the current generation of gravitational wave observations. This includes different NR approaches, modes and remnant properties. The distinction with the current models is also clear and those still show signs of systematic errors with respect to the most accurate NR waveforms.
2  2  0.9989  0.9990  0.9990  0.9347  244.54 
2  1  0.9965  0.9972  0.9968  0.6257  96.12 
2  0  0.9972  0.9973  0.9966  0.3091  56.06 
2  1  0.9982  0.9983  0.9983  0.5797  102.66 
2  2  0.9986  0.9986  0.9986  0.9600  215.48 
3  3  0.9901  0.9902  0.9912    29.63 
3  2  0.9887  0.9913  0.9902    17.14 
3  1  0.9785  0.9811  0.9801    8.98 
3  0  0.9803  0.9814  0.9834    5.57 
3  1  0.9848  0.9845  0.9847    9.17 
3  2  0.9867  0.9864  0.9862    17.01 
3  3  0.9899  0.9896  0.9901    28.87 
4  4  0.9921  0.9927  0.9938    11.99 
4  3  0.9800  0.9798  0.9814    6.61 
4  2  0.9830  0.9851  0.9838    4.17 
4  1  0.9856  0.9871  0.9868    2.30 
4  0  0.9317  0.9341  0.9377    1.55 
4  1  0.9854  0.9862  0.9861    2.32 
4  2  0.9825  0.9845  0.9836    4.26 
4  3  0.9827  0.9825  0.9835    6.93 
4  4  0.9906  0.9911  0.9919    10.13 
5  5  0.9703  0.9819  0.9848    3.19 
5  4  0.9646  0.9681  0.9735    1.72 
5  3  0.9641  0.9674  0.9708    1.09 
5  2  0.9575  0.9743  0.9765    0.66 
5  1  0.9657  0.9722  0.9734    0.36 
5  0  0.8730  0.8897  0.9013    0.25 
5  1  0.9636  0.9695  0.9710    0.37 
5  2  0.9541  0.9728  0.9765    0.67 
5  3  0.9688  0.9718  0.9738    1.13 
5  4  0.9643  0.9692  0.9725    1.71 
5  5  0.9657  0.9796  0.9825    2.73 
2  2  0.9993  282.02  0.9940  220.13  0.9993  259.04  0.9991  241.44  0.9996  237.10 

2  1  0.9985  28.16  0.9933  127.06  0.9973  96.60  0.9859  24.19  0.9991  19.24 
2  0  0.9294  5.81  0.9145  72.81  0.9975  54.88  0.9795  6.72  0.9816  5.14 
2  1  0.9986  28.13  0.9946  130.33  0.9985  103.84  0.9859  24.19  0.9991  19.24 
2  2  0.9993  282.02  0.9965  201.91  0.9990  226.98  0.9991  241.44  0.9996  237.10 
3  3  0.9463  6.49  0.9870  22.17  0.9904  31.88  0.9363  15.23  0.9978  37.33 
3  2  0.9993  7.80  0.9827  17.40  0.9915  17.89  0.9947  6.06  0.9986  6.15 
3  1  0.8844  1.17  0.9850  9.31  0.9788  9.30  0.7229  1.95  0.9099  0.73 
3  0  0.7757  0.78  0.8437  5.46  0.9792  5.73  0.8199  0.69  0.8720  0.28 
3  1  0.8879  1.21  0.9766  9.99  0.9838  9.51  0.7228  1.95  0.9099  0.73 
3  2  0.9993  7.80  0.9930  18.29  0.9861  17.81  0.9947  6.06  0.9986  6.15 
3  3  0.9463  6.49  0.9843  22.66  0.9898  30.67  0.9363  15.23  0.9978  37.33 
4  4  0.9922  12.17  0.9810  9.20  0.9926  12.97  0.9908  10.32  0.9956  13.75 
4  3  0.9930  2.24  0.9824  7.66  0.9794  7.11  0.9869  1.68  0.9967  1.94 
4  2  0.5968  0.69  0.9874  4.90  0.9841  4.37  0.9150  0.53  0.9068  0.49 
4  1  0.6361  0.14  0.9743  2.51  0.9857  2.39  0.7761  0.15  0.6647  0.10 
4  0  0.3514  0.49  0.6951  1.45  0.9355  1.60  0.2246  0.52  0.2578  0.28 
4  1  0.6533  0.13  0.9606  2.43  0.9848  2.45  0.7761  0.15  0.6647  0.10 
4  2  0.5967  0.69  0.9887  4.89  0.9830  4.51  0.9150  0.53  0.9068  0.49 
4  3  0.9930  2.24  0.9851  7.75  0.9830  7.37  0.9869  1.68  0.9967  1.94 
4  4  0.9922  12.17  0.9814  8.28  0.9912  10.78  0.9908  10.32  0.9956  13.75 
5  5  0.8662  0.64  0.9614  1.85  0.9709  3.29  0.9028  1.38  0.9807  3.83 
5  4  0.9758  0.58  0.9666  1.81  0.9648  1.85  0.9615  0.44  0.9893  0.61 
5  3  0.7287  0.12  0.9775  1.26  0.9631  1.16  0.7838  0.13  0.8748  0.19 
5  2  0.4258  0.22  0.9592  0.74  0.9552  0.71  0.6812  0.12  0.7639  0.13 
5  1  0.1107  0.07  0.9494  0.37  0.9658  0.38  0.6209  0.03  0.4149  0.04 
5  0  0.3735  0.10  0.5100  0.25  0.8761  0.27  0.3827  0.08  0.2603  0.06 
5  1  0.0553  0.04  0.9128  0.34  0.9639  0.40  0.6208  0.03  0.4152  0.04 
5  2  0.4258  0.22  0.9711  0.72  0.9509  0.74  0.6812  0.12  0.7639  0.13 
5  3  0.7286  0.12  0.9804  1.30  0.9683  1.22  0.7838  0.13  0.8748  0.19 
5  4  0.9758  0.58  0.9764  1.87  0.9661  1.86  0.9615  0.44  0.9893  0.61 
5  5  0.8662  0.64  0.9650  1.77  0.9663  2.76  0.9028  1.38  0.9807  3.83 
UID1N100  0.955294  0.619052  402.78 

UID1N118  0.955310  0.619079  404.40 
UID1N140  0.955311  0.619100  405.63 
UID1L2  0.955782  0.618905  
UID1L3  0.955813  0.618893  
UID1L4  0.955829  0.618899  
UID2N100  0.963445  0.581627  962.55 
UID2N118  0.963418  0.581480  996.49 
UID2N140  0.963405  0.581392  1016.06 
UID2L1  0.963768  0.579124  
UID2L2  0.964063  0.579988  
UID2L3  0.964063  0.579958  
UID3N100  0.961903  0.659634  614.70 
UID3N118  0.961920  0.659725  598.96 
UID3N140  0.961927  0.659781  587.63 
UID3L1  0.962123  0.658707  
UID3L2  0.962388  0.657731  
UID3L3  0.962401  0.657599  
UID4N100  0.962020  0.529128  312.81 
UID4N118  0.962028  0.529129  313.65 
UID4N140  0.962030  0.529130  314.05 
UID4L1  0.962114  0.528897  
UID4L2  0.962184  0.529023  
UID4L3  0.962174  0.529117  
UID5N100  0.968160  0.531761  171.57 
UID5N118  0.968171  0.531837  175.35 
UID5N140  0.968173  0.531873  177.81 
UID5L1  0.967872  0.531920  
UID5L2  0.968041  0.531934  
UID5L3  0.968051  0.531917 
UID  (RIT)  (SXS)  (GT)  (SEOB)  Model  

N100  N118  N140  L3  M120  (at NR)  
#1  60.4  61.0  61.0  60.9    62.7  v4 
#2  61.0  60.9  60.6  60.9    61.4  v3 
#3  60.4  60.5  60.7  60.7    60.4  v3 
#4  60.6  60.7  60.8  60.3  60.4  62.2  v4 
#5  60.0  60.0  60.1  60.0  59.8  61.2  v4 
NR  Label  Sim. ID  Model  

Group  (at NR)  
RIT  a  d0_D10.52_q1.3333_a0.25_n100 
0.7500  ( 0, 0, 0 )  ( 0, 0, 0.25 )  63.0  62.5  v4 
GT  b  (0.0,1.15) 
0.8696  ( 0, 0, 0 )  ( 0, 0, 0 )  62.2  61.5  v4 
RIT  c  q50_a0_a8_th_135_ph_30 
0.5000  ( 0, 0, 0 )  ( 0.490, 0.283, 0.566 )  62.5  60.7  v3 
BAM  d  BAM150914:24 
0.8912  ( 0.278, 0.605, 0.085 )  ( 0.151, 0.396, 0.017 )  62.7  61.0  v3 
SXS  e  SXS:BBH:0052 
0.3333  ( 0.001, 0.008, 0.499 )  ( 0.494, 0.073, 0.001 )  62.3  60.4  v3 
Iv Comparing NR simulations with observations of GW170104
The comparisons above demonstrate that our simulations agree with one another, but differ from approximate and phenomenological models oftused to describe precessing mergers. Fortunately, nature has provided us with a natural benchmark with which to assess the efficacy of our calculations and the significance of any discrepancies: observations of BHBH mergers. We use standard techniques Abbott et al. (The LIGO Scientific Collaboration and the Virgo Collaboration) (2016); Lange et al. (2017) to directly compare GW170104 to our simulations. For context, we also compare these observations to the corresponding predictions of approximate and phenomenological models that purport to describe the same event.
Figure 4 displays the direct comparison of the nonprecessing simulations by RIT and SXS complementary approaches for the configurations UID#1, UID#4, and UID#5, as given in Table 2. They directly compare to the signals as observed by LIGO H1 and L1 and with each other. The lower panel shows the residuals of the signals with respect to the RIT N118 simulation and compares it with the direct difference of the two approaches and also with the difference of the N118 and N100 resolutions, that measure the finite difference error of the N118, given the observed near 4th order convergence seen when included the N140 run into the analysis. For all three cases we note that the differences in any of the simulations is much smaller than the residuals and hence typical noise of the observations. This shows that the fast response runs performed to simulate BBH (lowmedium resolution) are in an acceptable good agreement with the expected higher resolution ones at the required level of errors.
Figure 5 displays the two precessing targeted simulations for GW170104 studied in this paper. We compare them with the L1 and H1 signals in grey and light grey in the plots. Here we also perform a double test of the accuracy of the simulations by considering the two main approaches to solve BBHs by the RIT and SXS groups and by considering the internal consistency of convergence of the waveforms with increasing resolutions. The waveforms again show a good agreement among themselves and their differences, shown in the lower panels are smaller than the residuals of the signals with respect to the N118 simulations. They are larger than in the aligned cases due to the choice of the initial spin configurations at at slightly different reference orbital frequencies. Note also that this comparisons do not align the peak of the waveforms and hence if independently fit to data would show much smaller differential residuals.
iv.1 Residuals versus resolution
For each UID, direct comparison of our simulations to the data selects a fiducial total mass which best fits the observations, as measured by the marginalized likelihood of the data assuming our simulations. Using the same mass for all simulations performed for that UID (e.g., by all groups and for all resolutions), we can for each simulation select the binary extrinsic parameters like event time and sky location which maximize the likelihood of the data, given our simulation and mass. Then, using these extrinsic parameters, we evaluate the expected detector response in the LIGO Hanford (H1) and Livingston (L1) instruments. This procedure has been used to reconstruct the gravitational wave signal for GW150914 Abbott et al. (The LIGO Scientific Collaboration and the Virgo Collaboration) (2016) and GW170104 events.
Figure 6 shows an example of these reconstructions for the highest LogLikelihood NR waveform (top candidate in Table 8 and UID#1. The top panel of this figure shows the NR predicted response in Hanford (bluered); the Hanford data (grey); and the Livingston data (in dark grey; shifted by 2.93ms and sign flipped). The bottom panel shows the residuals, and the difference between the two simulations in green. Note that the difference between waveforms is small compared to the residuals, but enough to make the simulation in blue (top candidate in Table 8) have a slightly higher Likelihood (63.0 vs. 62.5) over UID#1 in red, to match the signals over the whole range of frequencies considered. The same simulation resolution have been considered in both cases.
We have also analyzed the finite differences errors produced by fastresponse, low resolution (yet in the convergence regime) simulations of BBH mergers. The low, medium and high resolutions runs, N100, N118, and N140 respectively, by the RIT group show a nearly 4th order convergence (There are detailed studies of convergence for similar simulations in Refs. Healy et al. (2014); Healy and Lousto (2017); Healy et al. (2017c)) that allow to extrapolate to infinite resolution and evaluate the magnitude of the errors in the waveforms as compared to the residuals for this GW170104 event. We thus can evaluate the error of the N118 simulation is given by the (N100N118) difference, while the error of the N100 waveform is twice this difference and that of the N140 waveform is half that difference. This is displayed in the lower half of each panel in Figs. 4 and 5 and provide an alternative evaluation of the errors within a given NR method.
The studies carried out in this paper involving 3resolutions for each set of parameters well in the convergence regime of the simulations can be very costly from the resources point of view, totaling over 4 million service units (SUs) in computer clusters, as detailed in Table 9.
UID  N100  N118  N140  Total 

1  119  184  407  710 
2  313  451  557  1321 
3  145  217  476  838 
4  130  178  565  873 
5  67  118  306  491 
Total  774  1148  2311  4233 
According to the variations in the Table 7 that evaluates for the different resolutions we may derive as a rule of thumb that the N100 grid provides a good approximation for the nonprecessing binaries, while for the precessing ones, N118 is more appropriate. This leads to a reduction of the SUs needed for these 5 simulations, totaling nearly 1 million SUs, two thirds of which are due to the two precessing cases. The pseudospectral approach used by the SXS collaboration requires similar total wallclock times than the above finite differences approach, but spends an order of magnitude less resources. For instance, UID#1 (SXS:BBH:0626) required 11 kSUs for Lev4, 7.4 kSUs for Lev3, and 4.7 kSUs hours for Lev2.
iv.2 Likelihood of NR and models
For any proposed coalescing binary, characterized by its outgoing radiation as a function of all directions, we can compute a single quantity to assess its potential similarity to GW170104, accounting for all possible ways of orienting the source and placing it in the universe: the marginalized likelihood () Pankow et al. (2015); Abbott et al. (2016e); Lange et al. (2017). To provide a sense of scale, the distribution of over the posterior distribution including all intrinsic parameters is roughly universal Abbott et al. (2016e), approximately distributed as where has degrees of freedom (i.e., a mean value of , and its 90% credible interval is , where and for and , respectively). For each UID and for each proposed total mass , direct comparison of our simulations to the data allows us to compute a single number measuring the quality of fit: the marginalized likelihood . The maximum value of this function (here denoted by ) therefore measures the overall quality of fit. Table 7 shows for the five UIDs simulated here. For comparison, the last column shows calculated using an approximate model for the radiation from a coalescing binary. Obviously, if these approximate models and our simulations agree, then we should find the same result for at the same parameters. Finally, for context, the peak value of computed using SEOBNRv3 with generic parameters is 63.3. If our simulation parameters are wellchosen (and if both our simulations and these models are close to true solutions of Einstein’s equations), then this peak value should be in good agreement with the evaluated using our simulations.
First and foremost, up to Monte Carlo and fitting error, the marginalized likelihoods calculated with NR agree with each other comparing different resolutions and different approaches to solve the BBH problem, as required given the high degree of similarity between the underlying simulations. Second, the marginalized likelihoods computed at these proposed points are substantially below the largest found with approximate models like SEOBNRv3, except for UID3. Similar to the explanation described in Appendix B, the exception here is due to the differences between the precessing models ( was calculated with SEOB but the parameters were suggested with IMRP). Likewise, the binary parameters at which the peak value of occurs for SEOBNRv3 are substantially different from any of the proposed parameters explored here. This discrepancy suggests that the modelbased procedure that we adopted to target our followup simulations was not effective at finding the most likely parameters, as measured with . The poor performance of our targeted followup cannot simply reflect sampling error; even though the likelihood surface is nearly flat near the peak, so small errors are amplified in parameter space, this nearflatness also insures that systematic offsets should produce a small change in , if the underlying waveform calculations agree; see Appendix B for further discussion. Instead, we suspect the biases in arise because the models only approximate the correct solution of Einstein’s equations. Third, we confirm our hypothesis in Table 8 simply by demonstrating that other simulations (not performed in followup) fit the data substantially better than our targeted parameters.
On the one hand, NR followup simulations guided by the models (as displayed in Table 7) leads to lower marginalized likelihoods (). Conversely, other simulations shown in Table 8 produce higher , at points in parameter space where the models predict lower . This discrepancy suggest the two processes ( evaluated with NR and with the models) favor different regions of parameter space. In particular, table 8, which has one of the largest values of among all of the (roughly two thousand) simulations available to us, shows that the top precessing simulation is . This simulation has a mass ratio of 1:2, i.e. q=1/2, where the smaller hole is nonspinning and the larger hole is spinning with an intrinsic spin magnitude of 0.8 and pointing initially in a direction downwards with respect to the orbital angular momentum (=135 degrees) and an angle of 30 degrees from the line joining the two black holes (=30 degrees). This simulation belongs to a family of 6 simulations performed in Ref. Zlochower and Lousto (2015) labeled as NQ50TH135PH[0,30,60,90,120,150]. Those runs, supplemented by two control runs with angles we performed for this paper, are displayed in Fig. 7 versus the for this GW170104 event. The lower panels plots all those simulation with respect to their angle at merger as defined in Ref. Zlochower and Lousto (2015) and given in table XXI in that paper. The continuous curve provide a fit (detailed in table 10) for such values as reference and an estimate of the maximum value located near the phi=30 simulation.
0  0  62.3  54.9 
30  19.5  62.5  55.2 
60  34.8  62.2  54.1 
90  56.5  62.5  54.4 
120  98.5  61.6  54.1 
150  146.5  60.6  54.5 
210  194.7  59.3  55.1 
310  294.0  60.4  54.6 
0.38  
0.37 
The notable results displayed in Fig. 7, where seems to be sensitive to the orientation of the spin of the larger hole on the orbital plane, are consistent with broader trends that can be extracted using similar simulations: here, the set of 24 simulations of the family NQ50TH[30,60,90,135]PH[0,30,60,90,120,150] given in Ref. Zlochower and Lousto (2015) supplemented by the two aligned runs NQ50TH[0,180]PH0 given in Ref. Healy and Lousto (2017) and two runs specifically performed for this paper, NQ50TH135PH[200,310]. These simulations all have , a nonspinning smaller BH, and a spinning BH with fixed spin magnitude but changing orientation. Figure 8 shows a colormap derived from the maximum obtained for each of these simulations, using standard (MatLab) plotting tools. The last surface levels indicates the regions of largest likelihood (60,61,62) and a maximum, marked with an X, is located at TH=137, PH=87 with of 62.6. This results allow us to perform followup simulations seeking for this maximum.
In the plot, the black points are the NR simulations and the black curves are level sets of the colormap. Instead of plotting in the angles theta and phi, we plot in the HammerAitoff coordinates ^{2}^{2}2Weisstein, Eric W.“HammerAitoff EqualArea Projection.” From MathWorld–A Wolfram Web Resource. http://mathworld.wolfram.com/HammerAitoffEqualAreaProjection.html, which is a coordinate system where the whole angular space can be viewed as a 2d map. The points at the top left and bottom left are the poles, at the top, and at the bottom. The line connecting the two is the line. As you move from left to right, increases from 0 to 150 degrees (the maximum value of available in these simulations).
iv.3 Reconstructed NR waveforms
The analysis above – a difference in for models that should represent the same physical binary which is comparable to the expected range of over the posterior – suggests modest tension can exist between our NR simulations and the models used to draw inferences about GW170104. To illustrate this tension, in Fig. 9 we display the 90% confidence intervals of the precessing follow up cases (#2 and #3) computed by the two approximate/phenomenological models comparing them with the full numerical simulations (RIT’s with N100 resolutions, note that increasing the numerical resolutions to N118 and N140 reinforces this point). For each simulation, the waveform is generated by first fixing the total mass – selected by maximizing – and then choosing extrinsic parameters which maximize the likelihood. At merger, these reconstructed waveforms appear to be in modest tension with the confidence interval reported for ; for example, the peaks and troughs of the yellow (NR) curves are consistently at the boundaries of what the 90% credible intervals derived from waveforms allow. This illustration, however, relies on a non representative metric to assess waveform similarity (i.e, differences in the GW strain as a function of time, without reference to detector sensitivity, assessed by eye).
To remedy this deficiency, Figure 10 uses the match to compare our reconstructed NR waveforms with reconstructed waveforms drawn from the posterior parameter distribution of GW170104. The top panel uses a violin plot to illustrate the distribution of matches, with a solid bar showing the median value. The bottom panel shows a sample cumulative distribution. The median and maximum of this distribution provides a measure of how consistent the estimate via NR is with the distribution provided by the model. Using the maximum likelihood waveform from the model and posterior, these distributions should be proportional to a (centrally) distributed quantity, with median mismatch for the number of model degrees of freedom and the signal to noise ratio (SNR), where the specific choice for depends on the signal and detector/network being studied (e.g., for GW170104, the network SNR was ) . By contrast, in several of these overlap distributions, the peak and median values are manifestly offset downward, supporting a significant systematic difference between the radiation predicted from our approximate models and our NR waveforms, each generated from targeted NR followup simulations using parameters drawn from these selfsame model parameter distributions.
iv.4 Discussion
Using comparisons to data via as our guide, we found in Section IV.2 that modelbased and NRbased analyses seem to have maxima (in ) in different parts of parameter space; see Appendix B for greater detail. In the region identified as a good fit by modelbased analysis, corresponding NR simulations have a low . Conversely, several NR simulations with distinctly different parameters had a larger than the corresponding targeted NR simulations and modelbased comparisons evaluated at the same parameters. The two functions , evaluated using models and NR on parameters designed to be similar to and representative of plausible parameters for GW170104, do not agree, implying systematic differences between models and NR (i.e., a change in ). While we have for simplicity adopted one procedure which identifies candidate parameters to select our followup simulations, we emphasize that the specific procedure is largely arbitrary, as in this work we simply demonstrate the two marginalized likelihoods (NR and modelbased) disagree somewhere. Changes in the exact location and value of the marginalized likelihood are of less interest changes in the full posterior distribution; the latter subject is beyond the scope of this study.
The NR followup simulations and Bayesian inferences used in this work were performed soon after the identification of GW170104, and as such did not benefit from recent improvements in waveform modeling. Notably, by calibrating to a large suite of numerical relativity simulations, surrogate waveform models have been generated that, in a suitable part of parameter space, are markedly superior to any of the waveform models used for parameter inference to date Blackman et al. (2015, 2017b). Parameter inferences performed with these models should be more reliable and (by optimizing ) enable better targets for NR followup simulations.
For simplicity and brevity, we have directly compared our nonprecessing and precessing simulations to only one of the two extant families of phenomenological waveform models (SEOBv3/v4). While the two models are in good agreement for nonprecessing binaries, the other model (IMRP) has technical complications that limit its utility for our study. On the one hand, we cannot generate a similar waveform with similar initial conditions, preventing us from performing the straightforward comparisons shown in Figure 3. [As a frequency domain model, it did not adopt the same time conventions as NR and timedomain models for the precession phase (see,e.g., Williamson et al 2017 Williamson et al. (2017)).] On the other hand, the implementations available do not provide a spinweighted spherical harmonic decomposition, preventing us from performing the modebymode mismatch calculations in Table 4.
Previous investigations have demonstrated by example that posterior inferences with approximate waveform models can be biased, even for parameters consistent with observed binary black hole Williamson et al. (2017); Abbott et al. (2017). For example, a previous large study using simulations consistent with GW150914 found that, despite the brevity and relative simplicity of its signal, the inferred parameters could be biased for certain binary configurations relative to the line of sight Abbott et al. (2017), and much less so for others (e.g., nonprecessing and comparablemass binaries). The relevance and frequency of these configurations is not yet determined and depends on the binary black hole population which nature provides.
V Conclusions
After the detection of GW170104 Abbott et al. (2017c), we performed several simulations of binary black hole mergers, intending to reproduce LIGO’s observations using simulations with similar parameters. The parameters used were selected based on LIGO’s reported inferences about GW170104, generated by comparing two approximate models for binary black hole merger to the GW170104 data. Comparing these targeted simulations of binary black hole mergers, we find good agreement. We have shown that the differences among typical numerical simulations, used as a measure of their error, is much smaller (by over an order of magnitude) than the residuals of observation versus theory. On the other hand, we demonstrate (expected) differences between our numerical solutions to general relativity and the approximate models used to target our simulations. Because we used these models to identify candidate parameters for followup, our followup simulations were systematically biased away from the bestfitting parameters. These biases are not surprising, as the models used do not fully incorporate all the physics of binary merger, including higher modes and all features of precession, and are known to modestly disagree both with one another and with NR simulations. This does not mean that the models are not recovering the full signal: both models and NR could find similar likelihoods, but for different parameters. These bias can be particularly large for small mass ratios and highly spinning precessing binaries. We demonstrate that other, preexisting simulations with different parameters fit the data substantially better than the configurations targeted by modelbased techniques.
We have shown here (and in previous studies Bustillo et al. (2015); Lange et al. (2017)) that the standard low resolution, fastresponse, simulations provide an accurate description of GW signals, and can improve over the parameters determined by the models (See Table 8 and Fig. 7) for precessing and nonprecessing cases (note that while SEOBNRv4 improves on the inaccurate Lovelace et al. (2016) SEOBNRv2 Taracchini et al. (2014), it is still not at comparable accuracy to the NR simulations, See Figs. 23, for instance). The tension between the models and the full numerical simulations (notwithstanding Abbott et al. (2016f)) may be crucial in determining parameters such as individual spin of the holes and tests of general relativity for the large SNR signals, where the limitations of the models is larger). Both this study, focused on GW170104, and the investigation by Williamson et al. (2017), carried out on GW151226, point to the limitations of existing models to accurately determine binary parameters in the case of precessing BBH.
Regarding prospects for future followups, Figure 11 shows the distributions of the minimal total mass of the BBH systems in the NR catalogs Mroue et al. (2013); Jani et al. (2016); Healy et al. (2017a) given a starting gravitational wave frequency of 20 or 30 Hz in the source frame and its cumulative. This provides a coverage for the current events observed by LIGO (redshift effects improve this coverage by a factor of , where is the redshift). Coverage of lower total masses would require longer simulations or hybridization of the current waveforms.
Finally, we demonstrated the power of using purely numerical waveforms to determine parameters of a binary black hole merger as the previous case of GW150914 Abbott et al. (2016e) and similarly in the case of the source GW170104. More work is needed though to systematically and robustly include hybridization of waveforms and the case of generically precessing binaries Lange et al. (2017).
Vi Acknowledgments
The authors thank S. Husa, C. Berry, K.Chatziioannou, and A. Zimmerman for their feedback on this work. The RIT authors gratefully acknowledge the NSF for financial support from Grants PHY1607520, No. PHY1707946, No. ACI1550436, No. AST1516150, No. ACI1516125, No. PHY1726215. This work used the Extreme Science and Engineering Discovery Environment (XSEDE) [allocation TGPHY060027N], which is supported by NSF grant No. ACI1548562. Computational resources were also provided by the NewHorizons and BlueSky Clusters at the Rochester Institute of Technology, which were supported by NSF grants No. PHY0722703, No. DMS0820923, No. AST1028087, and No. PHY1229173. ROS is supported by NSF AST1412449, PHY1505629, and PHY1607520. The GT authors gratefully acknowledge NSF support through awards PHY1505824, 1505524, and XSEDE TGPHY120016. They also gratefully acknowledge support from CullenPeck and Dunn Families. The SXS collaboration authors gratefully acknowledge the NSF for financial support from Grants: No. PHY1307489, No. PHY1606522, PHY1606654, and AST 1333129. They also gratefully acknowledge support for this research at CITA from NSERC of Canada, the Ontario Early Researcher Awards Program, the Canada Research Chairs Program, and the Canadian Institute for Advanced Research. Calculations were done on the ORCA computer cluster, supported by NSF grant PHY1429873, the Research Corporation for Science Advancement, CSU Fullerton, the GPC supercomputer at the SciNet HPC Consortium Loken et al. (2010); SciNet is funded by: the Canada Foundation for Innovation (CFI) under the auspices of Compute Canada; the Government of Ontario; Ontario Research Fund (ORF) – Research Excellence; and the University of Toronto. Further calculations were performed on the Briarée cluster at Sherbrooke University, managed by Calcul Québec and Compute Canada and with operation funded by the Canada Foundation for Innovation (CFI), Ministére de l’Économie, de l’Innovation et des Exportations du Quebec (MEIE), RMGA and the Fonds de recherche du Québec  Nature et Technologies (FRQNT). PK gratefully acknowledges support from the Sherman Fairchild Foundation and NSF grants PHY1606654 and AST 1333129 at Cornell.
References
 Abbott et al. (2016a) B. P. Abbott et al. (LIGO Scientific Collaboration, Virgo Collaboration), Phys. Rev. Lett. 116, 061102 (2016a), arXiv:1602.03837 [grqc] .
 Abbott et al. (2016b) B. P. Abbott et al. (LIGO Scientific Collaboration, Virgo Collaboration), Phys. Rev. Lett. 116, 241103 (2016b), arXiv:1606.04855 [grqc] .
 Abbott et al. (2016c) B. P. Abbott et al. (Virgo, LIGO Scientific), Phys. Rev. X6, 041015 (2016c), arXiv:1606.04856 [grqc] .
 Abbott et al. (2017) B. P. Abbott, R. Abbott, T. D. Abbott, F. Acernese, K. Ackley, C. Adams, T. Adams, P. Addesso, R. X. Adhikari, V. B. Adya, and et al., Physical Review Letters 118, 221101 (2017), arXiv:1706.01812 [grqc] .
 Abbott et al. (2017a) B. P. Abbott et al. (Virgo, LIGO Scientific), (2017a), arXiv:1711.05578 [astroph.HE] .
 Abbott et al. (2017b) B. P. Abbott et al. (Virgo, LIGO Scientific), Phys. Rev. Lett. 119, 141101 (2017b), arXiv:1709.09660 [grqc] .
 Taracchini et al. (2014) A. Taracchini, A. Buonanno, Y. Pan, T. Hinderer, M. Boyle, D. A. Hemberger, L. E. Kidder, G. Lovelace, A. H. Mroué, H. P. Pfeiffer, M. A. Scheel, B. Szilágyi, N. W. Taylor, and A. Zenginoglu, Phys. Rev. D 89, 061502 (2014), arXiv:1311.2544 [grqc] .
 Pürrer (2014) M. Pürrer, Class. Quantum Grav. 31, 195010 (2014), arXiv:1402.4146 [grqc] .
 Hannam et al. (2014) M. Hannam, P. Schmidt, A. Bohé, L. Haegel, S. Husa, et al., Phys. Rev. Lett. 113, 151101 (2014), arXiv:1308.3271 [grqc] .
 Abbott et al. (The LIGO Scientific Collaboration and the Virgo Collaboration) (2016) B. Abbott et al. (The LIGO Scientific Collaboration and the Virgo Collaboration), Phys. Rev. D 94, 064035 (2016).
 Lange et al. (2017) J. Lange, R. O’Shaughnessy, M. Boyle, J. Calderón Bustillo, M. Campanelli, T. Chu, J. A. Clark, N. Demos, H. Fong, J. Healy, D. Hemberger, I. Hinder, K. Jani, B. Khamesra, L. E. Kidder, P. Kumar, P. Laguna, C. O. Lousto, G. Lovelace, S. Ossokine, H. Pfeiffer, M. A. Scheel, D. Shoemaker, B. Szilagyi, S. Teukolsky, and Y. Zlochower, ArXiv eprints (2017), arXiv:1705.09833 [grqc] .
 Pretorius (2005) F. Pretorius, Phys. Rev. Lett. 95, 121101 (2005), arXiv:grqc/0507014 [grqc] .
 Campanelli et al. (2006) M. Campanelli, C. O. Lousto, P. Marronetti, and Y. Zlochower, Phys. Rev. Lett. 96, 111101 (2006), arXiv:grqc/0511048 [grqc] .
 Baker et al. (2006) J. G. Baker, J. Centrella, D.I. Choi, M. Koppitz, and J. van Meter, Phys. Rev. Lett. 96, 111102 (2006), arXiv:grqc/0511103 [grqc] .
 Mroue et al. (2013) A. H. Mroue et al., Phys. Rev. Lett. 111, 241104 (2013), arXiv:1304.6077 [grqc] .
 Jani et al. (2016) K. Jani, J. Healy, J. A. Clark, L. London, P. Laguna, and D. Shoemaker, Class. Quant. Grav. 33, 204001 (2016), arXiv:1605.03204 [grqc] .
 Healy et al. (2017a) J. Healy, C. O. Lousto, Y. Zlochower, and M. Campanelli, (2017a), arXiv:1703.03423 [grqc] .
 Lange et al. (2017) J. Lange et al., (2017), arXiv:1705.09833 [grqc] .
 Healy et al. (2014) J. Healy, C. O. Lousto, and Y. Zlochower, Phys. Rev. D 89, 104052 (2014), arXiv:1406.7295 [grqc] .
 Ghosh et al. (2016) A. Ghosh et al., Phys. Rev. D94, 021101 (2016), arXiv:1602.02453 [grqc] .
 Abbott et al. (2016d) B. P. Abbott et al. (LIGO Scientific Collaboration, Virgo Collaboration), Phys. Rev. Lett. 116, 221101 (2016d), arXiv:1602.03841 [grqc] .
 Zlochower et al. (2005) Y. Zlochower, J. Baker, M. Campanelli, and C. Lousto, Phys. Rev. D 72, 024021 (2005), arXiv:grqc/0505055 [grqc] .
 Marronetti et al. (2008) P. Marronetti, W. Tichy, B. Bruegmann, J. Gonzalez, and U. Sperhake, Phys. Rev. D 77, 064010 (2008), arXiv:0709.2160 [grqc] .
 Lousto and Zlochower (2008) C. O. Lousto and Y. Zlochower, Phys. Rev. D 77, 024034 (2008), arXiv:0711.1165 [grqc] .
 Löffler et al. (2012) F. Löffler, J. Faber, E. Bentivegna, T. Bode, P. Diener, R. Haas, I. Hinder, B. C. Mundim, C. D. Ott, E. Schnetter, G. Allen, M. Campanelli, and P. Laguna, Class. Quantum Grav. 29, 115001 (2012), arXiv:1111.3344 [grqc] .
 (26) Einstein Toolkit home page: http://einsteintoolkit.org.
 (27) Cactus Computational Toolkit home page: http://cactuscode.org.
 Schnetter et al. (2004) E. Schnetter, S. H. Hawley, and I. Hawke, Class. Quant. Grav. 21, 1465 (2004), grqc/0310042 .
 Healy et al. (2017b) J. Healy, C. O. Lousto, H. Nakano, and Y. Zlochower, Class. Quant. Grav. 34, 145011 (2017b), arXiv:1702.00872 [grqc] .
 Brandt and Brügmann (1997) S. Brandt and B. Brügmann, Phys. Rev. Lett. 78, 3606 (1997), grqc/9703066 .
 Ansorg et al. (2004) M. Ansorg, B. Bruegmann, and W. Tichy, Phys. Rev. D70, 064011 (2004), arXiv:grqc/0404056 [grqc] .
 Thornburg (2004) J. Thornburg, Class. Quant. Grav. 21, 743 (2004), grqc/0306056 .
 Dreyer et al. (2003) O. Dreyer, B. Krishnan, D. Shoemaker, and E. Schnetter, Phys. Rev. D67, 024018 (2003), grqc/0206008 .
 Campanelli et al. (2007) M. Campanelli, C. O. Lousto, Y. Zlochower, B. Krishnan, and D. Merritt, Phys. Rev. D 75, 064030 (2007), arXiv:grqc/0612076 [grqc] .
 Campanelli and Lousto (1999) M. Campanelli and C. O. Lousto, Phys. Rev. D 59, 124022 (1999), arXiv:grqc/9811019 [grqc] .
 Lousto and Zlochower (2007) C. O. Lousto and Y. Zlochower, Phys. Rev. D 76, 041502 (2007), arXiv:grqc/0703061 [GRQC] .
 Nakano et al. (2015) H. Nakano, J. Healy, C. O. Lousto, and Y. Zlochower, Phys. Rev. D 91, 104022 (2015), arXiv:1503.00718 [grqc] .
 Healy et al. (2017c) J. Healy, C. O. Lousto, and Y. Zlochower, Phys. Rev. D96, 024031 (2017c), arXiv:1705.07034 [grqc] .
 Brügmann et al. (2008) B. Brügmann, J. A. González, M. Hannam, S. Husa, U. Sperhake, and W. Tichy, Phys. Rev. D 77, 024027 (2008), grqc/0610128 .
 Herrmann et al. (2007a) F. Herrmann, I. Hinder, D. M. Shoemaker, P. Laguna, and R. A. Matzner, Phys. Rev. D76, 084032 (2007a), arXiv:0706.2541 [grqc] .
 Herrmann et al. (2007b) F. Herrmann, I. Hinder, D. Shoemaker, P. Laguna, and R. A. Matzner, Astrophys. J. 661, 430 (2007b), arXiv:grqc/0701143 [GRQC] .
 Hinder et al. (2008) I. Hinder, B. Vaishnav, F. Herrmann, D. Shoemaker, and P. Laguna, Phys. Rev. D77, 081502 (2008), arXiv:0710.5167 [grqc] .
 Healy et al. (2009a) J. Healy, F. Herrmann, I. Hinder, D. M. Shoemaker, P. Laguna, and R. A. Matzner, Phys. Rev. Lett. 102, 041101 (2009a), arXiv:0807.3292 [grqc] .
 Hinder et al. (2010) I. Hinder, F. Herrmann, P. Laguna, and D. Shoemaker, Phys. Rev. D 82, 024033 (2010), arXiv:0806.1037 [grqc] .
 Healy et al. (2009b) J. Healy, J. Levin, and D. Shoemaker, Phys. Rev. Lett. 103, 131101 (2009b), arXiv:0907.0671 [grqc] .
 Healy et al. (2010) J. Healy, P. Laguna, R. A. Matzner, and D. M. Shoemaker, Phys. Rev. D 81, 081501 (2010), arXiv:0905.3914 [grqc] .
 Bode et al. (2010) T. Bode, R. Haas, T. Bogdanovic, P. Laguna, and D. Shoemaker, Astrophys. J. 715, 1117 (2010), arXiv:0912.0087 [grqc] .
 Husa et al. (2006) S. Husa, I. Hinder, and C. Lechner, Comput. Phys. Commun. 174, 983 (2006), arXiv:grqc/0404023 [grqc] .
 (49) http://www.blackholes.org/SpEC.html.
 Lovelace et al. (2008) G. Lovelace, R. Owen, H. P. Pfeiffer, and T. Chu, Phys. Rev. D 78, 084017 (2008).
 York (1999) J. W. York, Phys. Rev. Lett. 82, 1350 (1999).
 Pfeiffer et al. (2003) H. P. Pfeiffer, L. E. Kidder, M. A. Scheel, and S. A. Teukolsky, Comput. Phys. Commun. 152, 253 (2003), grqc/0202096 .
 Ossokine et al. (2015) S. Ossokine, F. Foucart, H. P. Pfeiffer, M. Boyle, and B. Szilágyi, Class. Quantum Grav. 32, 245010 (2015), arXiv:1506.01689 [grqc] .
 Caudill et al. (2006) M. Caudill, G. B. Cook, J. D. Grigsby, and H. P. Pfeiffer, Phys. Rev. D 74, 064011 (2006), grqc/0605053 .
 Pfeiffer et al. (2007) H. P. Pfeiffer, D. A. Brown, L. E. Kidder, L. Lindblom, G. Lovelace, and M. A. Scheel, Class. Quantum Grav. 24, S59 (2007), grqc/0702106 .
 Buonanno et al. (2011) A. Buonanno, L. E. Kidder, A. H. Mroué, H. P. Pfeiffer, and A. Taracchini, Phys. Rev. D 83, 104034 (2011), arXiv:1012.1549 [grqc] .
 Mroué and Pfeiffer (2012) A. H. Mroué and H. P. Pfeiffer, (2012), arXiv:1210.2958 [grqc] .
 Lindblom et al. (2006) L. Lindblom, M. A. Scheel, L. E. Kidder, R. Owen, and O. Rinne, Class. Quantum Grav. 23, S447 (2006), grqc/0512093 .
 Friedrich (1985) H. Friedrich, Commun. Math. Phys. 100, 525 (1985).
 Garfinkle (2002) D. Garfinkle, Phys. Rev. D 65, 044029 (2002).
 Pretorius (2005) F. Pretorius, Class. Quantum Grav. 22, 425 (2005), grqc/0407110 .
 Lindblom and Szilágyi (2009) L. Lindblom and B. Szilágyi, Phys. Rev. D 80, 084019 (2009), arXiv:0904.4873 .
 Choptuik and Pretorius (2010) M. W. Choptuik and F. Pretorius, Phys. Rev. Lett. 104, 111101 (2010), arXiv:0908.1780 [grqc] .
 Szilágyi et al. (2009) B. Szilágyi, L. Lindblom, and M. A. Scheel, Phys. Rev. D 80, 124010 (2009), arXiv:0909.3557 [grqc] .
 Lovelace et al. (2011) G. Lovelace, M. A. Scheel, and B. Szilágyi, Phys. Rev. D 83, 024010 (2011), arXiv:1010.2777 [grqc] .
 Szilágyi (2014) B. Szilágyi, Int. J. Mod. Phys. D 23, 1430014 (2014), arXiv:1405.3693 [grqc] .
 M. A. Scheel, M. Boyle, T. Chu, L. E. Kidder, K. D. Matthews and H. P. Pfeiffer (2009) M. A. Scheel, M. Boyle, T. Chu, L. E. Kidder, K. D. Matthews and H. P. Pfeiffer, Phys. Rev. D 79, 024003 (2009), arXiv:grqc/0810.1767 .
 Hemberger et al. (2013) D. A. Hemberger, M. A. Scheel, L. E. Kidder, B. Szilágyi, G. Lovelace, N. W. Taylor, and S. A. Teukolsky, Class. Quantum Grav. 30, 115001 (2013), arXiv:1211.6079 [grqc] .
 Ossokine et al. (2013) S. Ossokine, L. E. Kidder, and H. P. Pfeiffer, Phys. Rev. D 88, 084031 (2013), arXiv:1304.3067 [grqc] .
 Scheel et al. (2015) M. A. Scheel, M. Giesler, D. A. Hemberger, G. Lovelace, K. Kuper, M. Boyle, B. Szilágyi, and L. E. Kidder, Class. Quantum Grav. 32, 105009 (2015), arXiv:1412.1803 [grqc] .
 Rinne (2006) O. Rinne, Class. Quantum Grav. 23, 6275 (2006).
 Rinne et al. (2007) O. Rinne, L. Lindblom, and M. A. Scheel, Class. Quantum Grav. 24, 4053 (2007).
 Gundlach (1998) C. Gundlach, Phys. Rev. D 57, 863 (1998).
 Cook and Whiting (2007) G. B. Cook and B. F. Whiting, Phys. Rev. D 76, 041501(R) (2007).
 Owen (2007) R. Owen, Topics in Numerical Relativity: The periodic standingwave approximation, the stability of constraints in free evolution, and the spin of dynamical black holes, Ph.D. thesis, California Institute of Technology (2007).
 Boyle et al. (2007) M. Boyle, D. A. Brown, L. E. Kidder, A. H. Mroué, H. P. Pfeiffer, M. A. Scheel, G. B. Cook, and S. A. Teukolsky, Phys. Rev. D 76, 124038 (2007), arXiv:0710.0158 [grqc] .
 Buchman and Sarbach (2007) L. T. Buchman and O. C. A. Sarbach, Class. Quantum Grav. 24, S307 (2007).
 Rinne et al. (2009) O. Rinne, L. T. Buchman, M. A. Scheel, and H. P. Pfeiffer, Class. Quantum Grav. 26, 075009 (2009).
 Bishop et al. (1997) N. T. Bishop, R. Gomez, L. Lehner, M. Maharaj, and J. Winicour, Phys. Rev. D 56, 6298 (1997), arXiv:grqc/9708065 .
 Winicour (2005) J. Winicour, Living Rev. Rel. 8 (2005).
 Gomez et al. (2007) R. Gomez, W. Barreto, and S. Frittelli, Phys. Rev. D 76, 124029 (2007), arXiv:0711.0564 [grqc] .
 Reisswig et al. (2013) C. Reisswig, N. T. Bishop, and D. Pollney, Gen. Rel. Grav. 45, 1069 (2013), arXiv:1208.3891 [grqc] .
 Handmer and Szilágyi (2015) C. J. Handmer and B. Szilágyi, Class. Quantum Grav. 32, 025008 (2015), arXiv:1406.7029 .
 Boyle and Mroué (2009) M. Boyle and A. H. Mroué, Phys. Rev. D 80, 124045 (2009), arXiv:0905.3177 [grqc] .
 Taylor et al. (2013) N. W. Taylor, M. Boyle, C. Reisswig, M. A. Scheel, T. Chu, L. E. Kidder, and B. Szilágyi, Phys. Rev. D 88, 124010 (2013), arXiv:1309.3605 [grqc] .
 Chu et al. (2016) T. Chu, H. Fong, P. Kumar, H. P. Pfeiffer, M. Boyle, D. A. Hemberger, L. E. Kidder, M. A. Scheel, and B. Szilágyi, Class. Quantum Grav. 33, 165001 (2016), arXiv:1512.06800 [grqc] .
 Boyle (2016a) M. Boyle, Phys. Rev. D 93, 084031 (2016a).
 Abbott et al. (2017c) B. P. Abbott et al. (VIRGO, LIGO Scientific), Phys. Rev. Lett. 118, 221101 (2017c), arXiv:1706.01812 [grqc] .
 Veitch et al. (2015) J. Veitch, V. Raymond, B. Farr, W. Farr, P. Graff, S. Vitale, B. Aylott, K. Blackburn, N. Christensen, M. Coughlin, W. Del Pozzo, F. Feroz, J. Gair, C.J. Haster, V. Kalogera, T. Littenberg, I. Mandel, R. O’Shaughnessy, M. Pitkin, C. Rodriguez, C. Röver, T. Sidery, R. Smith, M. Van Der Sluys, A. Vecchio, W. Vousden, and L. Wade, Phys. Rev. D 91, 042003 (2015).

(90)
Https://dcc.ligo.org/DocDB/0123/T1500606/006/
NRInjectionInfrastructure.pdf.  Boyle (2016b) M. Boyle, Phys. Rev. D93, 084031 (2016b), arXiv:1509.00862 [grqc] .
 Kelly and Baker (2013) B. J. Kelly and J. G. Baker, Phys. Rev. D87, 084004 (2013), arXiv:1212.5553 [grqc] .
 Lovelace et al. (2016) G. Lovelace, C. O. Lousto, J. Healy, M. A. Scheel, A. Garcia, R. O’Shaughnessy, M. Boyle, M. Campanelli, D. A. Hemberger, L. E. Kidder, H. P. Pfeiffer, B. Szilágyi, S. A. Teukolsky, and Y. Zlochower, Submitted to Class. Quantum Grav.; arXiv:1607.05377 (2016), arXiv:1607.05377 [grqc] .
 Bohé et al. (2017) A. Bohé, L. Shao, A. Taracchini, A. Buonanno, S. Babak, I. W. Harry, I. Hinder, S. Ossokine, M. Pürrer, V. Raymond, T. Chu, H. Fong, P. Kumar, H. P. Pfeiffer, M. Boyle, D. A. Hemberger, L. E. Kidder, G. Lovelace, M. A. Scheel, and B. Szilágyi, Phys. Rev. D 95, 044028 (2017), arXiv:1611.03703 [grqc] .
 Healy and Lousto (2017) J. Healy and C. O. Lousto, Phys. Rev. D95, 024037 (2017), arXiv:1610.09713 [grqc] .
 JimńezForteza et al. (2017) X. JimńezForteza, D. Keitel, S. Husa, M. Hannam, S. Khan, and M. Pürrer, Phys. Rev. D95, 064024 (2017), arXiv:1611.00332 [grqc] .
 Blackman et al. (2017a) J. Blackman, S. E. Field, M. A. Scheel, C. R. Galley, C. D. Ott, M. Boyle, L. E. Kidder, H. P. Pfeiffer, and B. Szilágyi, Phys. Rev. D96, 024058 (2017a), arXiv:1705.07089 [grqc] .
 London et al. (2017) L. London, S. Khan, E. FauchonJones, X. J. Forteza, M. Hannam, S. Husa, C. Kalaghatgi, F. Ohme, and F. Pannarale, (2017), arXiv:1708.00404 [grqc] .
 Pankow et al. (2015) C. Pankow, P. Brady, E. Ochsner, and R. O’Shaughnessy, Phys. Rev. D92, 023002 (2015), arXiv:1502.04370 [grqc] .
 Abbott et al. (2016e) B. P. Abbott et al. (LIGO Scientific Collaboration, Virgo Collaboration), Phys. Rev. D 94, 064035 (2016e), arXiv:1606.01262 [grqc] .
 Zlochower and Lousto (2015) Y. Zlochower and C. O. Lousto, Phys. Rev. D92, 024022 (2015), arXiv:1503.07536 [grqc] .
 (102) Weisstein, Eric W.“HammerAitoff EqualArea Projection.” From MathWorld–A Wolfram Web Resource. http://mathworld.wolfram.com/HammerAitoffEqualAreaProjection.html.
 Blackman et al. (2015) J. Blackman, S. E. Field, C. R. Galley, B. Szilágyi, M. A. Scheel, M. Tiglio, and D. A. Hemberger, Phys. Rev. Lett. 115, 121102 (2015), arXiv:1502.07758 [grqc] .
 Blackman et al. (2017b) J. Blackman, S. E. Field, M. A. Scheel, C. R. Galley, D. A. Hemberger, P. Schmidt, and R. Smith, (2017b), arXiv:1701.00550 [grqc] .
 Williamson et al. (2017) A. R. Williamson, J. Lange, R. O’Shaughnessy, J. A. Clark, P. Kumar, J. Calderón Bustillo, and J. Veitch, (2017), arXiv:1709.03095 [grqc] .
 Abbott et al. (2017) B. P. Abbott, R. Abbott, T. D. Abbott, M. R. Abernathy, F. Acernese, K. Ackley, C. Adams, T. Adams, P. Addesso, R. X. Adhikari, and et al., Classical and Quantum Gravity 34, 104002 (2017), arXiv:1611.07531 [grqc] .
 Bustillo et al. (2015) J. C. Bustillo, A. Bohé, S. Husa, A. M. Sintes, M. Hannam, et al., (2015), arXiv:1501.00918 [grqc] .
 Taracchini et al. (2014) A. Taracchini, A. Buonanno, Y. Pan, T. Hinderer, M. Boyle, D. A. Hemberger, L. E. Kidder, G. Lovelace, A. H. Mroue, H. P. Pfeiffer, M. A. Scheel, B. Szilágyi, N. W. Taylor, and A. Zenginoglu, Phys. Rev. D 89 (R), 061502 (2014), arXiv:1311.2544 [grqc] .
 Abbott et al. (2016f) B. P. Abbott et al. (Virgo, LIGO Scientific), (2016f), arXiv:1611.07531 [grqc] .
 Loken et al. (2010) C. Loken, D. Gruner, L. Groer, R. Peltier, N. Bunn, M. Craig, T. Henriques, J. Dempsey, C.H. Yu, J. Chen, L. J. Dursi, J. Chong, S. Northrup, J. Pinto, N. Knecht, and R. V. Zon, J. Phys.: Conf. Ser. 256, 012026 (2010).
Appendix A Impact of discrete posterior sampling
In the text, we performed several calculations that depend on on an inferred posterior distribution: identifying parameters for NR followup calculations; maximizing the marginalized likelihood ; and calculating mismatches. These calculations are performed using a finitesize collection of approximately independent, identicallydistributed samples from a posterior distribution Veitch et al. (2015). In this appendix we briefly quantify the (small) effects our finite sample size has on our conclusions and comparisons. For simplicity, we will conservatively standardize our calculations to posterior samples; in practice, usually many more were used.
The match and marginalized likelihood distributions are welldescribed by a distribution with a suitable number of degrees of freedom, corresponding to the model dimension of the intrinsic parameter space (i.e., for calculations which omit precession, and for calculations which include it). For example if is the true maximum marginalized likelihood, then the marginalized likelihood distribution over the posterior is wellapproximated by the distribution of where is distributed with degrees of freedom. If we have independent draws from the distribution, the smallest value of will be distributed according to , where is the cumulative distribution for the distribution. At 95% confidence, the maximum value of over the samples is therefore . As a result, if we estimate the maximum value of with the maximum over our posterior samples, we find am estimate which is smaller than the true maximum value by . Evaluating this expression for and in the conservative limit of only samples, we find a systematic sampling error of () in (), respectively, in our estimate of the peak marginalized likelihood. This systematic sampling error in our estimate of the peak marginalized likelihood is smaller than the differences in marginalized likelihoods discussed in the text and figures.
Likewise, given the number of samples, the targeted parameters should be very close to the true maximum a posteriori values. Qualitatively speaking, due to finite sample size effects, our estimate of each parameter has an uncertainty of roughly , or roughly of the width of the distribution using our fiducial sample size.
Appendix B Mixed messages: Maximum likelihood, , and a posteriori
One goal of this work is to demonstrate, by a concrete counterexample, that NR followup must be targeted and assessed selfconsistently.
One source of inconsistency in our original NR followup strategy was the algorithm by which NR simulations were selected from modelbased inference. Our NR followup simulations were selected by (approximately) maximizing the a posteriori probability, proportional to the (15dimensional) likelihood ; the (7dimensional) prior for extrinsic parameters ; and the (8dimensional) prior for intrinsic parameters . This maximum a posteriori (MaP) location does not generally correspond to the parameters which maximize the likelihood (maxL). The intrinsic parameters selected by both approaches also do not cause the marginalized likelihood to take on its largest value. In principle, to avoid introducing artificial inconsistencies simply due to the choice of point estimate, we should have targeted followup simulations using . To assess how much our choice of targeting impacted our estimate of , we evaluated the marginalized likelihood at our estimates of all three locations. Each location was approximated by our (finitesize) set of posterior samples. For the posterior distribution adopted to generate UID4 – a nonprecessing productionquality analysis where SEOBNRv4 was both used to generate the reference posterior used to find the MaP and maxL parameters and to compute a modelbased – we find that the modelbased values at the MaP and maxL points to be effectively indistinguishable due to Monte Carlo error (61.4 and 61.2 respectively, with an estimated Monte Carlo error of ). This similarity suggests that, when a fully selfconsistent analysis is performed, then even if the MaP and maxL parameters differ slightly, they will produce similar values of , with differences far smaller than the differences between NR and modelbased analysis.
For the reasons described in Section IV.4, we consistently adopt SEOBbased models to evaluate our modelbased . Because different phenomenological approximants do not agree, the posterior distributions used to identify the parameters for UID3 and 5 used an IMRD/Pbased approximant produce different MaP and maxL parameters. Conversely, to the extent these models agree, they should estimate model parameters corresponding to the same values of . In fact, however, when we evaluate with SEOBNRv4 on the MaP and maxL parameters of the posterior used to find UID5, we find both values disagree with the values seen for UID4, being lower (60.8) and higher (62) respectively. These differences in clearly indicate small differences between the two modelbased analyses, comparable to (but smaller than) the differences seen between modelbased analyses and NR.