Targeted numerical simulations of binary black holes for GW170104

Targeted numerical simulations of binary black holes for GW170104

J. Healy    J. Lange    R. O’Shaughnessy    C.O. Lousto    M. Campanelli    A. R. Williamson    Y. Zlochower Center for Computational Relativity and Gravitation, School of Mathematical Sciences, Rochester Institute of Technology, 85 Lomb Memorial Drive, Rochester, New York 14623    J. Calderón Bustillo    J.A. Clark    C. Evans    D. Ferguson    S. Ghonge    K. Jani    B  Khamesra    P. Laguna    D. M. Shoemaker Center for Relativistic Astrophysics and School of Physics, Georgia Institute of Technology, Atlanta, GA 30332, USA    M. Boyle Center for Astrophysics and Planetary Science, Cornell University, Ithaca, New York 14853, USA    A. García Gravitational Wave Physics and Astronomy Center, California State University Fullerton, Fullerton, California 92834, USA    D. A. Hemberger Theoretical Astrophysics 350-17, California Institute of Technology, Pasadena, CA 91125, USA    L. E. Kidder Center for Astrophysics and Planetary Science, Cornell University, Ithaca, New York 14853, USA    P. Kumar Center for Astrophysics and Planetary Science, Cornell University, Ithaca, New York 14853, USA Canadian Institute for Theoretical Astrophysics, 60 St. George Street, University of Toronto, Toronto, ON M5S 3H8, Canada    G. Lovelace Gravitational Wave Physics and Astronomy Center, California State University Fullerton, Fullerton, California 92834, USA    H. P. Pfeiffer Max Planck Institute for Gravitational Physics (Albert Einstein Institute), Am Mühlenberg 1, 14476 Potsdam-Golm, Germany Canadian Institute for Theoretical Astrophysics, 60 St. George Street, University of Toronto, Toronto, ON M5S 3H8, Canada    M. A. Scheel Theoretical Astrophysics 350-17, California Institute of Technology, Pasadena, CA 91125, USA    S.A. Teukolsky Center for Astrophysics and Planetary Science, Cornell University, Ithaca, New York 14853, USA Theoretical Astrophysics 350-17, California Institute of Technology, Pasadena, CA 91125, USA

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 LIGO-Virgo 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 state-of-the-art 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 best-fitting 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 merger-ringdown signals as tests of general relativity Ghosh et al. (2016); Abbott et al. (2016d, c).

ii.1 Simulations using finite-difference, moving-puncture 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, sixth-order finite differencing in space Lousto and Zlochower (2008) and a fourth-order Runge Kutta time integrator (the code does not upwind the advection terms) and a 5th-order Kreiss-Oliger 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 post-Newtonian 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#1-5, 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. Sixth-order 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 Kerr-Schild 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 quasi-equilibrium 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 first-order 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 adaptively-refined Lovelace et al. (2011); Szilágyi (2014) spatial grid that extends from pure-outflow 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 time-stepping 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. Constraint-preserving 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 fast-flow 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 Newman-Penrose 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 Regge-Wheeler-Zerilli-Moncrief equations at finite radius Buchman and Sarbach (2007); Rinne et al. (2009), and Cauchy-Characteristic 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 Regge-Wheeler-Zerilli-Moncrief 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 well-motivated 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
Table 1: Follow-up Parameter Table

Spin Conventions: are specified in a frame where (i) , i.e. the Newtonian orbital angular momentum is along the z-axis. (ii) the vector pointing from to is the x-axis, (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 spin-components given below are those consistent with Eqs. (9)-(11) of 111

The label, UID#1-5, 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 spin-vectors 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 NR-frequency 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 set-up 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 after-junk masses/spins, ideally extracted at the reference frequency. For precessing runs, in particular, the spin-components 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

#1 -7.9168 6.7407 -2.829e-4 0.07467 0.3196 0.3056 0.4599 0.5401 0.9928 0.7343 -0.8267 6e-4
#2 -7.8211 4.1029 -4.837e-4 0.07837 0.3277 0.4400 0.3441 0.6559 0.9922 0.1977 0.7580 1e-3
#3 -8.7720 4.2543 -3.316e-4 0.07160 0.2796 0.3584 0.3266 0.6734 0.9930 0.5101 0.8445 1e-3
#4 -8.4742 6.0567 -2.918e-4 0.07421 0.3991 0.4219 0.4168 0.5832 0.9928 0.2205 -0.7110 2e-4
#5 -9.4395 4.0593 -2.718e-4 0.06698 0.2753 0.6865 0.3007 0.6993 0.9933 -0.3634 -0.1256 5e-4
Table 2: Initial data parameters for the quasi-circular configurations with a smaller mass black hole (labeled 1), and a larger mass spinning black hole (labeled 2). The punctures are located at and , with momenta , mass parameters , horizon (Christodoulou) masses , total ADM mass , dimensionless spins , and eccentricity, .
UID1-N100 0.459757 0.539780 0.851718 (0,0,0.737493) (0,0,-0.829064) 0.002586
UID1-N118 0.459758 0.539801 0.851718 (0,0,0.737473) (0,0,-0.829040) 0.002587
UID1-N140 0.459758 0.539801 0.851718 (0,0,0.737464) (0,0,-0.829030) 0.002587
UID1-L2 0.459913 0.540183 0.851401 (0,0,0.734254) (0,0,-0.827442) 0.003135
UID1-L3 0.459902 0.540171 0.851400 (0,0,0.734138) (0,0,-0.827657) 0.003139
UID1-L4 0.459907 0.540176 0.851402 (0,0,0.734106) (0,0,-0.827690) 0.003139
UID2-N100 0.344090 0.655693 0.524773 ( 0.119494, 0.144676, -0.063539 ) ( -0.675724, 0.181910, -0.293427 ) 0.004010
UID2-N118 0.344091 0.655693 0.524774 ( 0.118992, 0.144981, -0.063537 ) ( -0.676089, 0.180798, -0.293141 ) 0.004010
UID2-N140 0.344091 0.655694 0.524774 ( 0.118754, 0.145124, -0.063541 ) ( -0.676262, 0.180273, -0.293005 ) 0.004010
UID2-L1 0.344082 0.655876 0.524614 ( 0.101171, 0.156396, -0.066135 ) ( -0.693687, 0.096203, -0.290389 ) 0.004001
UID2-L2 0.344069 0.655966 0.524523 ( 0.100390, 0.156863, -0.066508 ) ( -0.693803, 0.095200, -0.289578 ) 0.004013
UID2-L3 0.344072 0.655958 0.524534 ( 0.100596, 0.156736, -0.066591 ) ( -0.693727, 0.095442, -0.289439 ) 0.004010
UID3-N100 0.326605 0.672963 0.485325 ( 0.360515, -0.132028, -0.338700 ) ( -0.685165, 0.491281, -0.068777 ) 0.003308
UID3-N118 0.326606 0.672963 0.485326 ( 0.360458, -0.131982, -0.338637 ) ( -0.685154, 0.491237, -0.068765 ) 0.003308
UID3-N140 0.326607 0.672964 0.485326 ( 0.360429, -0.132025, -0.338591 ) ( -0.685106, 0.491276, -0.068763 ) 0.003308
UID3-L1 0.326576 0.673382 0.484978 ( 0.358221, -0.095241, -0.350438 ) ( -0.749729, 0.387741, -0.052141 ) 0.003330
UID3-L2 0.326598 0.673451 0.484962 ( 0.358434, -0.094398, -0.350478 ) ( -0.750204, 0.386567, -0.051977 ) 0.003339
UID3-L3 0.326583 0.673462 0.484931 ( 0.352293, -0.114550, -0.350800 ) ( -0.727233, 0.427282, -0.059549 ) 0.003308
UID4-N100 0.416817 0.583036 0.714908 (0,0,0.221608) (0,0,-0.712189) 0.002629
UID4-N118 0.416819 0.583037 0.714911 (0,0,0.221602) (0,0,-0.712184) 0.002630
UID4-N140 0.416820 0.583037 0.714912 (0,0,0.221600) (0,0,-0.712182) 0.002630
UID4-L1 0.416790 0.583195 0.714666 (0,0,0.220660) (0,0,-0.710930) 0.002635
UID4-L2 0.416809 0.583214 0.714675 (0,0,0.220453) (0,0,-0.710967) 0.002635
UID4-L3 0.416817 0.583203 0.714703 (0,0,0.220427) (0,0,-0.710930) 0.002619
UID5-N100 0.300721 0.699282 0.430043 (0,0,-0.366067) (0,0,-0.125360) 0.002932
UID5-N118 0.300722 0.699284 0.430043 (0,0,-0.366031) (0,0,-0.125352) 0.002932
UID5-N140 0.300723 0.699285 0.430043 (0,0,-0.366015) (0,0,-0.125348) 0.002932
UID5-L1 0.300697 0.699261 0.430021 (0,0,-0.363371) (0,0,-0.125641) 0.002900
UID5-L2 0.300708 0.699223 0.430061 (0,0,-0.363424) (0,0,-0.125550) 0.002864
UID5-L3 0.300722 0.699178 0.430107 (0,0,-0.363345) (0,0,-0.125616) 0.002855
Table 3: Values of the individual masses, mass ratio, dimensionless spins, and frequency, given at a time after the gauge settles. For the nonprecessing cases (1, 4, and 5), this time is for RIT and for SXS. For the precessing systems (2 and 3), the values are given such that the relaxed frequency is the same between RIT and SXS.

For followup #1, the initial spurious burst contains a non-negligible 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. (4-5) in Boyle (2016b) and Eqs. (7-8) in Kelly and Baker (2013). Note that this does not change the physical waveforms, only how they are spread over modes.

Figure 1: The small BH (top) and large BH (bottom) spins for followup #3. RIT’s simulation has solid lines; SXS’s has dashed lines; and spin evolution as predicted by SEOBNRv3 is shown with a wide-dashed line. The spin components, , , and , are red, blue, and green respectively.

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 Bowen-York 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 observationally-driven measure: the overlap or match. The black and grey lines in Figures 2 and 3 show the match between the two simulations’ (RIT-SXS and RIT-GT 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 time-shift is computed by maximizing the overlap of the (2,2) mode; this phase- and time-shift 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 mode-by-mode, for all resolutions.

Based on Figures 2, 3 and Table 4 we anticipate the lowest production-quality 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 mode-by-mode 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 over-resolved at these resolutions. Those remnant quantities are also important to model fitting formulae Healy and Lousto (2017); Jimńez-Forteza 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
Table 4: Match between individual spherical harmonic modes of the SXS and RIT UID3 waveforms, using the H1 PSD characterizing data near GW170104. Following Lovelace et al. (2016), rather than maximize over time and phase for each independently, our mode-by-mode comparisons fix the event time and overall phase using one mode (here, the (2,2) mode). The successively higher resolution simulations from RIT, labeled as are compared to the L3 (highest) resolution run from SXS. The minimal frequency is taken as for and for for a fiducial total mass of . The column labeled shows the match between RIT N140 and the corresponding SEOBNRv3 mode. Rows with a “-” are not modeled by SEOBNRv3. The column labeled shows the overlap of N140 with itself, with a minimum frequency of in all cases, to indicate the significance of the mode.
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
Table 5: Matches () between individual spherical harmonic modes of the SXS and RIT waveforms, using the H1 PSD characterizing data near GW170104. The lowest resolution simulations from RIT, labeled N100, are compared to the L3 resolution run from SXS. The minimal frequency is taken as for and for . The columns labeled show the overlap of N100 with itself, , to indicate the significance of the mode.
UID1-N100 0.955294 0.619052 402.78
UID1-N118 0.955310 0.619079 404.40
UID1-N140 0.955311 0.619100 405.63
UID1-L2 0.955782 0.618905
UID1-L3 0.955813 0.618893
UID1-L4 0.955829 0.618899
UID2-N100 0.963445 0.581627 962.55
UID2-N118 0.963418 0.581480 996.49
UID2-N140 0.963405 0.581392 1016.06
UID2-L1 0.963768 0.579124
UID2-L2 0.964063 0.579988
UID2-L3 0.964063 0.579958
UID3-N100 0.961903 0.659634 614.70
UID3-N118 0.961920 0.659725 598.96
UID3-N140 0.961927 0.659781 587.63
UID3-L1 0.962123 0.658707
UID3-L2 0.962388 0.657731
UID3-L3 0.962401 0.657599
UID4-N100 0.962020 0.529128 312.81
UID4-N118 0.962028 0.529129 313.65
UID4-N140 0.962030 0.529130 314.05
UID4-L1 0.962114 0.528897
UID4-L2 0.962184 0.529023
UID4-L3 0.962174 0.529117
UID5-N100 0.968160 0.531761 171.57
UID5-N118 0.968171 0.531837 175.35
UID5-N140 0.968173 0.531873 177.81
UID5-L1 0.967872 0.531920
UID5-L2 0.968041 0.531934
UID5-L3 0.968051 0.531917
Table 6: Remnant results for spinning binaries. We show the remnant mass in units of the total initial mass , the remnant dimensionless spin , and the remnant velocity in the x-y plane . We show results for different LazEv resolutions (N100, N118, and N140) and different SpEC resolutions (L0, L2, L4, and L6).
Figure 2: For the three nonprecessing UIDs # 1,4,5 in Table 2, matches between SXS, RIT, and SEOBNRv4 (2,2) modes as a function of , using the H1 PSD characterizing data near GW170104. We also compare with GT runs for UIDs # 4,5. Compare to also to similar plots for GW150914 Lovelace et al. (2016).
Figure 3: For the two precessing UIDs#2,3 in Table 2, matches between SXS, RIT, and SEOBNRv3 (2,2) modes as a function of as a function of , using the H1 PSD characterizing data near GW170104. In this comparison, the mode of all three simulations and SEOBNRv3 are extracted relative to the axis, identified from their common initial orbital parameters. While these frame identifications are coordinate-dependent for precessing binaries – implying our comparisons here could include both intrinsic disagreement and systematic error due to (say) overall misalignment – the good agreement shown in Figure 1 for the equally coordinate-dependent spins suggests that convention-dependent sources contribute little to the mismatches illustrated here.
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
Table 7: Marginalized likelihood of the data: This table shows the results for the 5 simulations when directly compared to the data. For these results, we use the same PSD adopted in all other calculations, with Hz (i.e. low-frequency cutoff). The first column is the UID. The second column is the estimated peak log marginalized likelihood , maximized over binary total mass, for the NR followup simulation. The third column is the corresponding log marginalized likelihood, using exactly the same intrinsic parameters (e.g., masses and spins) as maximize the likelihood in the second column, evaluated using a phenomenological approximate model instead of numerical relativity. The fourth column is the specific model used: either SEOBNRv3 (for precessing simulations) or SEOBNRv4 (for nonprecessing simulations). To see more on this parameter estimation method, see Abbott et al. (The LIGO Scientific Collaboration and the Virgo Collaboration) (2016); Lange et al. (2017).
NR Label Sim. ID Model
Group (at NR)
RIT a d0_D10.52_q1.3333_a-0.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
Table 8: Marginalized likelihood of the data: Selected other simulations: This table shows the results for several other simulations that particularly match the data well and the SEOB model results at those parameter points. These simulations are part of the top 15 simulations in . When comparing the NR values here to the ones in Table 7, one can see these to be generally higher i.e. better match the data. When comparing the NR values to the SEOB at the same points, one sees a consistent lower SEOB value This implies that these points were not picked for NR Followup due to the lower SEOB value.

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 oft-used 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 BH-BH 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 (low-medium resolution) are in an acceptable good agreement with the expected higher resolution ones at the required level of errors.

Figure 4: Comparison of the GW170104 signal seen by LIGO detectors H1 and L1 (in grey and dark grey) with the computer simulations of black hole mergers from SXS, RIT, and GT approaches for the nonprecessing cases labeled as , , and in Table 2.

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.

Figure 5: Comparison of the GW170104 signal seen by LIGO detectors H1 and L1 (in grey and dark grey) with the computer simulations of black hole mergers from SXS and RIT approaches for the nonprecessing cases labeled as and in Table 2.

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 Log-Likelihood NR waveform (top candidate in Table 8 and UID#1. The top panel of this figure shows the NR predicted response in Hanford (blue-red); 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 fast-response, 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 (N100-N118) 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 3-resolutions 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
Table 9: kSUs (1000 core-hours) for each RIT run and resolution.

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 well-chosen (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 model-based 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 near-flatness 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.

Figure 6: Comparison of the GW170104 signal seen by LIGO detectors H1 and L1 (in grey and dark grey) with the computer simulations of black hole mergers from RIT at low resolution for the nonprecessing case labeled as in Table 2 and the highest value for an NR simulation given in Table 8 (d0_D10.52_q1.3333_a-0.25_n100).

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.

Figure 7: The log-likelihood of the NQ50TH135 series Zlochower and Lousto (2015) assuming a period of versus initial angle (top panel) and merger angle (bottom panel.) Data (red) and fits (blue) are given in Table 10.
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
Table 10: The log-likelihood of the NQ50TH135 series Zlochower and Lousto (2015). Fittings of the form is also given for both the initial and .

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 color-map 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.

Figure 8: The log-likelihood of the NQ50THPHI series Zlochower and Lousto (2015) as a color map with red giving the highest and blue the lowest. The black dots (and grey diamonds, obtained by symmetry) represent the NR simulations and we have used Hammer-Aitoff coordinates , to represent the map and level curves with the top values of . The maximum, marked with an X, is located at TH=137, PH=87 reaching .

In the plot, the black points are the NR simulations and the black curves are level sets of the color-map. Instead of plotting in the angles theta and phi, we plot in the Hammer-Aitoff coordinates 222Weisstein, Eric W.“Hammer-Aitoff Equal-Area Projection.” From MathWorld–A Wolfram Web Resource., 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

Figure 9: Comparison of the 90% confidence intervals of GW170104 from the two precessing models with the computer simulations of black hole mergers (in orange) from the best-fitting NR simulations listed in Table 8.

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).

Figure 10: Distribution of the overall match between each NR waveform (a,b,c,d,e) listed in Table 8 and the (distribution of) waveforms produced by model-based parameter inference, as reported in Abbott et al. (2017c). Matches are produced for the signal in H1, L1, and overall. If an NR signal is perfectly consistent with these models for some parameters, then the distribution of matches will be well-approximated by a distribution with degrees of freedom. Many of the best-fitting NR simulations have a distribution of matches that is significantly offset relative to this expected distribution, reflecting the mild tension shown in Figure 9.

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 model-based and NR-based 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 model-based 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 model-based 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 model-based) 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 time-domain 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 spin-weighted spherical harmonic decomposition, preventing us from performing the mode-by-mode 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 comparable-mass 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 best-fitting 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, pre-existing simulations with different parameters fit the data substantially better than the configurations targeted by model-based techniques.

We have shown here (and in previous studies Bustillo et al. (2015); Lange et al. (2017)) that the standard low resolution, fast-response, 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 non-precessing 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. 2-3, 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.

Figure 11: This plot shows histograms of the of simulations in a given total mass bin of size , for starting at 20Hz and for starting at 30Hz with the given total mass. This is for all the runs in the public SXS+GT+RIT Catalogs Mroue et al. (2013); Jani et al. (2016); Healy et al. (2017a). The bottom histograms shows how many simulations in the catalogs can be used from 20Hz or 30Hz from a minimal mass on, i.e. the cumulative of the upper plot.

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 PHY-1607520, No. PHY-1707946, No. ACI-1550436, No. AST-1516150, No. ACI-1516125, No. PHY-1726215. This work used the Extreme Science and Engineering Discovery Environment (XSEDE) [allocation TG-PHY060027N], which is supported by NSF grant No. ACI-1548562. Computational resources were also provided by the NewHorizons and BlueSky Clusters at the Rochester Institute of Technology, which were supported by NSF grants No. PHY-0722703, No. DMS-0820923, No. AST-1028087, and No. PHY-1229173. ROS is supported by NSF AST-1412449, PHY-1505629, and PHY-1607520. The GT authors gratefully acknowledge NSF support through awards PHY-1505824, 1505524, and XSEDE TG-PHY120016. They also gratefully acknowledge support from Cullen-Peck and Dunn Families. The SXS collaboration authors gratefully acknowledge the NSF for financial support from Grants: No. PHY-1307489, No. PHY-1606522, PHY-1606654, 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 PHY-1429873, 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 (FRQ-NT). PK gratefully acknowledges support from the Sherman Fairchild Foundation and NSF grants PHY-1606654 and AST- 1333129 at Cornell.


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 finite-size collection of approximately independent, identically-distributed 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 well-described 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 well-approximated 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 self-consistently.

One source of inconsistency in our original NR followup strategy was the algorithm by which NR simulations were selected from model-based inference. Our NR followup simulations were selected by (approximately) maximizing the a posteriori probability, proportional to the (15-dimensional) likelihood ; the (7-dimensional) prior for extrinsic parameters ; and the (8-dimensional) 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 (finite-size) set of posterior samples. For the posterior distribution adopted to generate UID4 – a nonprecessing production-quality analysis where SEOBNRv4 was both used to generate the reference posterior used to find the MaP and maxL parameters and to compute a model-based – we find that the model-based 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 self-consistent 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 model-based analysis.

For the reasons described in Section IV.4, we consistently adopt SEOB-based models to evaluate our model-based . Because different phenomenological approximants do not agree, the posterior distributions used to identify the parameters for UID3 and 5 used an IMRD/P-based 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 model-based analyses, comparable to (but smaller than) the differences seen between model-based analyses and NR.

Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description