The discovery of GW150914_095045 by the Advanced LIGO [1] and Virgo [2] detectors marked the beginning of gravitational-wave (GW) astronomy. Since then, during the first three observing runs (referred to as O1, O2, and O3) [3][4][5][6], the LVK Collaboration, consisting of LIGO and Virgo, together with the KAGRA detector [7], has detected a total of 90 gravitational wave events [6]. Several independent analyses of the public data have further revealed another ∼ 15 events [8][9][10][11]. Analyses of the detected signals allow us to infer individual and population properties of merging compact objects [3][4][5][8][9][10][11][12][13][14][15][16][17][18][19][20] including the existence of mass gaps and possible formation channels of binary black hole (BBH) systems [21][22][23][24]. GW signals also encode information of the equation of state of the neutron-star matter [25][26][27][28], the nature of gravity [29][30][31][32][33], environments around compact objects [34], and inform our understanding of cosmology by providing an independent measurement of the Hubble constant [35][36][37].
Inference of binary source properties from a detected signal relies on the availability of expedient and accurate gravitational waveform models that are capable of describing the entire coalescence from the inspiral through the binary merger and ringdown of the remnant object. Current GW models can be categorized into three sets: phenomenological models [38][39][40][41][42][43][44][45][46][47][48][49], effective-one-body (EOB) models [50][51][52][53][54][55][56][57][58][59][60][61][62][63], and numerical relativity (NR) surrogate models [64][65][66][67][68][69]. The phenomenological and EOB waveform families are semi-analytical models that rely on a physically motivated ansatz, with calibration parameters that are fit to a set of NR simulations. NR surrogate mod-els instead take a purely data-driven approach by training the model directly on NR simulations without the need for an ansatz. NR surrogates have been shown to be more accurate than semi-analytical models within their common regime of validity and capture the underlying physics in the NR simulations (such as precession) at an accuracy comparable to the simulations themselves [66][67][68]. Consequently, NR surrogate models can provide highly accurate and trustworthy information about BBH source properties, especially as detector sensitivity improves [70][71][72][73].
Compared to phenomenological and EOB waveform models, however, surrogate models have a restricted regime of validity as they rely on the availability of NR simulations for training. The surrogate model we consider in this work, NRSur7dq4 [67], can only be evaluated for mass ratios q ≥ 1/6, where q := m 2,det /m 1,det ; m 1,det and m 2,det are the masses of the primary and secondary black hole, respectively, in the detector frame. Furthermore, because of the computational limitations of NR simulations, NRSur7dq4 only supports relatively short duration signals (∼ 20 orbits), which further restricts its applicability to heavy binaries with a detector-frame total mass of M det := m 1,det + m 2,det ≥ 60M ⊙ These restrictions leave 47 events from GWTC-3 that can be analyzed with the NRSur7dq4 waveform model (see Section III for details).
In this paper, we present a thorough analysis of these 47 events using the NRSur7dq4 waveform model for binary source parameter inference and an associated NRSur7dq4Remnant model to infer the mass, spin vector, and kick vector of the final black hole remnant. Our work builds upon previous studies that applied NR surrogate models to special events such as the first observation [70], the first observation of a high-mass ratio BBH merger [71], a highly precessing system [74,75], and an intermediatemass black hole formed through a BBH merger [73]. We also compare (i) NRSur7dq4 posteriors against publicly available LVK posteriors [5,6,76,77] obtained using the IMRPhenomXPHM [44] and SEOBNRv4PHM [55] models, (ii) the recovered Bayes factors and differences in the signal-to-noise ratios (SNRs) between NRSur7dq4 and IMRPhenomXPHM, and (iii) the remnant mass and spin magnitude estimates obtained for the three models. We find that while, for most events, our results are consistent with the LVK analyses, more than 20% of events exhibit noticeably different measurements and some differences, such as constraining GW191109_010717's effective spin to be confidently negative, may have important astrophysical implications.
The rest of the paper is organized as follows. Section II presents our analysis framework, including strain data, the NRSur7dq4 model, and our choice of priors, reference frame, and sampler settings. In Sec. III, we discuss how the set of 47 events that are analyzed have been selected. We provide an overview of the NRSur7dq4 results in Sec. IV and quantify the difference between NRSur7dq4 and IMRPhenomXPHM/SEOBNRv4PHM posteriors. The events that show the largest discrepancy with the LVK results are considered in more detail in Sec. IV A and Sec. IV C. In Sec. V, we use Bayes factors and recovered SNRs to better understand whether the data has a preference for a particular model, finding the data shows mild support for the NRSur7dq4 model. Constraints on the mass, spin vector, and kick vector are considered in Sec. VI. Finally, in Sec. VII we summarize our results and discuss the implications of our findings. Our posterior samples are publicly available at Ref. [78].
In this section, we summarize the Bayesian inference methods used in this study (Sec. II A) and provide an overview of our data analysis framework. Our setup is mostly consistent with settings used to obtain the public LVK posteriors [5,6,76,77], with the relevant differences detailed in Sec. II E and II F
The GW signal from a coalescing quasi-circular compact BBH system in general relativity can be completely characterized by a set of 15 parameters that we denote Θ = {I, E}. The parameter vector Θ consists of eight intrinsic parameters I and seven extrinsic parameters E. The vector I := {m 1 , m 2 , χ 1 , χ 2 , θ 1 , θ 2 , ϕ 12 , ϕ JL } contains the intrinsic parameters that describe the binary: the component masses m 1 and m 2 (with m 1 ≥ m 2 ), dimensionless spin magnitudes χ 1 and χ 2 , spin tilt angles θ 1 and θ 2 , and two azimuthl spin angles ϕ 12 and ϕ JL . The definition of these parameters is nicely summarized in Appendix E of Ref. [79]. We further use χ 1 and χ 2 to denote spin vectors for each component black hole, L to denote the orbital angular momentum, and J to denote the total angular momentum. Moreover, we distinguish between the source-frame and detector-frame masses by using the subscript 'det'. For example, m 1 is the source-frame mass and m 1,det the detector-frame mass of the primary black hole. These masses are related by m 1 = m 1,det /(1 + z) where z is the redshift of the source. The set of extrinsic parameters E = {θ JN , D L , α, δ, ψ, ϕ c , t c } parameterizes the location and orientation of the binary relative to the detectors as well as the coalescence time. The angle between the total angular momentum of the binary and the line-of-sight to the detector is denoted by θ JN (filling the role of the 'inclination angle') while the luminosity distance to the source is denoted by D L . Right ascension α and declination δ parameterize the location of the source in the sky, and ψ is the polarization angle. The coalescence time is denoted by t c while ϕ c indicates a reference orbital phase.
We use Bayes' theorem to compute the posterior probability distribution of the binary parameters,
where H is the signal hypothesis, d k (t) represents timedomain strain data for the k th detector,
which is assumed to be a sum of the true signal h k (t; Θ) and noise n k (t) in each detector. The subscript "k" on our signal model denotes that the observed signal will look different in different interferometers. The posterior p(Θ|{d k }, H) is the target for the parameter estimation analysis while the model evidence,
is the target for hypothesis testing (sometimes referred to as model selection). The prior probability, π(Θ|H), is a prescribed probability distribution and represents our initial assumptions about the parameters describing an individual binary. The likelihood function,
) describes how well each set of Θ matches the data. Here, h H k (Θ) is the signal waveform generated from a specific waveform model as part of our hypothesis (we typically omit the superscript for brevity), and ⟨a|b⟩ is the noiseweighted overlap integral defined as
with S n (f ) being the one-sided power spectral density (PSD) of the detector noise, a "∼" indicates a Fourier transform operation, and * represents the complex conjugate. The integration limits, f low and f high , are chosen to reflect the sensitivity bandwidth of the detectors; specific values are given in Sec. II F.
To quantify how much more likely that the data is described by a signal and not a noise process, we compute the Bayes factor,
where Z H = Z({d k }|H) and Z n denote, respectively, the evidence for a signal model H and a noise-only model.
Surrogate models use NR waveforms as training data and build a highly accurate interpolant over the parameter space using a combination of reduced-order modeling [80,81], parametric fits, and non-linear transformations of the waveform data [65]. The NRSur7dq4 model [67] used in this paper is trained on 1528 NR simulations [67] and spans the 7-dimensional parameter space of spin-precessing binaries. The model includes all ℓ ≤ 4 spin-weighted spherical harmonic modes, as well as the precession frame dynamics and spin evolution of the black holes. While the model has been trained for mass ratio 1/4 ≤ q ≤ 1.0 and spins 0.0 ≤ χ 1,2 ≤ 0.8, it can be extrapolated to q ≥ 1/6 and χ 1,2 ≤ 1.0 [67,82]. In its regime of validity, NRSur7dq4 improves upon semi-analytical models by about an order of magnitude in accuracy, in terms of mismatches against NR waveforms [67,82].
We compare our results with public LVK posteriors [5,6] obtained using the IMRPhenomXPHM [44,83,84] and SEOBNRv4PHM [55,85] waveform models.
IMRPhenomXPHM is a phenomenological model that includes the effects of precession and the {(ℓ, m)} = {(2, ±2), (2, ±1), (3, ±3), (3, ±2), (4, ±4)} modes in the coprecessing frame. SEOBNRv4PHM is an EOB model that also includes the effects of precession and the {(ℓ, m)} = {(2, ±2), (2, ±1), (3, ±3), (4, ±4), (5, ±5)} modes in the coprecessing frame modes. Importantly, while both IMRPhenomXPHM and SEOBNRv4PHM are calibrated against nonprecessing NR simulations, they are not calibrated against precessing simulations. Instead, they apply a frame-twisting procedure to approximately transform an aligned-spin waveform in the coprecessing frame to a precessing waveform in the inertial frame.
To infer the mass, spin vectors, and the kick velocity vector of the remnant black hole, we use the remnant surrogate model, NRSur7dq4Remnant [67,86]. The NRSur7dq4Remnant model is trained on the same set of 1528 NR simulations as NRSur7dq4, and is applicable in the same region of parameter space. NRSur7dq4Remnant improves upon previous remnant models by an order of magnitude in accuracy, in terms of errors in remnant properties with respect to NR simulations [67].
All time-dependent binary parameters, such as the BH spins and the system's orientation, are defined in a frame such that at some reference time (t ref ), the z-axis is along the instantaneous angular momentum vector L, the xaxis is along the line of separation from the less massive BH to the more massive BH, and the y-axis completes the right-handed triad. The reference point is defined as the time (t ref ) during the binary evolution where the GW frequency in the coprecessing frame (defined as twice the time derivative of Eq. 3 of Ref. [67]) passes f ref = 20 Hz. Following Ref. [87], we refer to this frame as the wave frame at f ref = 20 Hz. We measure all binary parameters at this reference frequency.
The strain d(t), PSD, and detector calibration uncertainty data are obtained from the LVK public data releases [76,77,[88][89][90]. We use the de-glitched strain data Events with SEOBNRv4PHM Events with de-glitched posteriors missing in strain data public LVK datasets
Table I. List of events for which SEOBNRv4PHM posteriors are absent in the publicly released LVK posteriors [5,6,76,77] as well as the events for which we use publicly available deglitched strain data [89,90,95].
for certain events as summarized in Table I. The publicly available PSDs were generated using the on-source BayesWave method [91][92][93] while the effect of frequencydependent uncertainties in amplitude and phase of the interferometer calibration on the parameter estimation of each event follows the methods of Refs. [79,94]. The signal's geocentric trigger time, duration, and sampling rate are also taken from publicly available 1 data releases. In particular, the event trigger times that can be separately obtained from the Gravitational Wave Open Science Center (GWOSC) were not used as they can be inconsistent 2 with the times used for the LVK parameter estimation analyses. The GWTC-3 data release [5,6,76,77,89,90] further describe the methods used for data conditioning and (where needed) deglitching.
Our assumptions for the priors π(Θ|H) are identical to the LVK analyses of the GWTC-3 catalog [3][4][5][6] with additional restrictions on the mass ratio and the total mass of the binary:
• We choose uniform priors in the detector-frame component masses subject to the following constraints 1 From each event's C01:IMRPhenomXPHM and C01:SEOBNRv4PHM parameter estimation results found within the respective PEDataRelease_mixed_cosmo.h5 files. 2 The geocentric trigger times in GWOSC come from a search pipeline, which we found to be not sufficiently precise given our time-of-coalescence priors, especially for multi-modal distributions in the time-of-arrival parameter.
Dynesty [97] live-points = 2000 tolerance = 0.1 nact = 50
Table II. Dynesty static nested sampling configuration parameters. We use 2000 live points and sample away from a current live point with random walks. We require that the number of random walk steps taken in each chain is at least 50 (=nact) times the autocorrelation length of the chain. Sampling continues until the estimated contribution of the remaining prior volume to the evidence is less than 0.1 (=tolerance).
on the BBH system: (i) the chirp mass satisfies 12M ⊙ ≤ M c,det ≤ 400M ⊙ , (ii) the mass ratio satisfies 0.167 ≤ q ≤ 1, and (iii) the total mass satisfies M det ≥ 60M ⊙ , where the chirp mass [96] is defined as
(m 1,det +m 1,det ) 1/5 . The second and third constraints are imposed to restrict the analysis to NRSur7dq4's region of validity; see Sec. III.
• Uniform priors are used for the dimensionless spin magnitudes (0.0 ≤ χ 1,2 ≤ 0.99) of the binary, with spin orientations taken as uniform on the unit sphere.
• The prior on the luminosity distance assumes uniform source distribution in comoving volume and time as implemented in the UniformSourceFrame prior class in Ref [79] within 100 Mpc ≤ D L ≤ 10, 000 Mpc. For some events, we use a higher upper bound to ensure that the posterior is contained within the prior's range.
• For the orbital inclination angle θ JN , we assume a uniform prior over -1 ≤ cos(θ JN ) ≤ 1.
• Priors on the sky location parameters, right ascension α and declination δ, are assumed to be uniform over the sky with periodic boundary conditions for α.
• For the time of coalescence, we assume a uniform prior over t geo -0.1 ≤ t c ≤ t geo + 0.1 where t geo is the geocentric trigger time.
Performing a parameter estimation analysis requires specifying numerous configuration options whose specific settings can impact the final inferred source parameters. We have chosen our analysis settings to match the LVK analyses of GWTC-3 as closely as possible based on what is documented in Ref. [88] and in the public LVK data release [76,77,89,90]. We briefly describe the most important parameter estimation settings while pointing out differences when they arise. Our complete configuration and environment files are made publicly available [78].
To estimate source properties, we employ the publicly-available Bayesian inference libraries parallel-bilby [98] (version 1.0.1) and bilby [79,99] (version 1.1.5) together with the dynesty [97] (version 1.0.1) nested sampling algorithm [100]. We report some of the most important sampler configuration settings in Table II. We have varied these values for several events to check our posteriors are sufficiently converged. We have also varied the number of processes used in our parallel-bilby computation. As documented in App. A 1, we notice that an important quantity to consider is live-points-per-process. When this number becomes too small, the computed posteriors are demonstrably inaccurate. Through extensive experimentation, we find that using 128 processes (corresponding to ≈ 16 live-points-per-process) provides robust posterior computations. Finally, following the official LVK analyses [76,77,89,90], for each event, we perform a total of four independent Bayesian inference runs with different initial random seeds to account for statistical randomness. We then combine the posteriors from these four analyses weighted by their individual evidences to compute the final posteriors. As an additional consistency check, for some events we have also compared posteriors with different initial seeds, finding no significant differences between these runs.
When computing the overlap integral in Eq.( 5), we follow Ref. [76,77,89,90] and set the maximum frequency to be f high = 0.875 × f PSD , where f PSD is the largest frequency in the publicly available PSD data [76,77]. The factor of 0.875 is used to minimize roll-off effects caused by the application of a tapering window to the time-domain data as implemented in bilby. We set the minimum frequency to be f low = 20Hz for most of the events, with notable exceptions being GW190521_030229 (f low = 11Hz for all detectors; following Ref. [101]) and GW190727_060333 (f low = 50Hz for the L1 detector; following Ref. [6] to exclude data that could be corrupted by the presence of a nearby glitch).
We call the NRSur7dq4 model through its LALSimulation interface [102].
For generating NRSur7dq4 waveforms, and as discussed in Sec. II C, we always use a reference frequency of f ref = 20Hz when setting model parameters. Furthermore, we generate the full length of the surrogate waveform (about 20 orbits) by passing f min = 0 to the waveform generator function. The important distinction between f low and f min is that f low sets the lower limit for the overlap integral (5), while f min sets the starting frequency of the waveform. By setting f min = 0 we obtain the longest possible signal from NRSur7dq4.
One important limitation of NRSur7dq4 is the restricted duration of the waveforms provided by the model. The model is only able to generate relatively short duration waveforms corresponding to about 20 orbits before merger, making it difficult to analyze low-mass systems (M det ≲ 60M ⊙ ) when using a lower cut-off frequency of f low = 20 Hz, the current default value for most LVK analyses [5,6,76,77]. Furthermore, the model is valid only for binaries with a mass ratio in the range 1/6 ≤ q ≤ 1. This restricts the number of events we can analyze with the NRSur7dq4 model.
We inspect the public LVK posteriors of the binary source properties obtained using IMRPhenomXPHM [44] and SEOBNRv4PHM [55] waveform models for all events in the GWTC-3 [5,6,76,77] catalog. We select the events for which at least 95% of the samples in the posteriors (inferred using the IMRPhenomXPHM model), have M det larger than 60M ⊙ and mass ratio q larger than 0.167. This ensures that these events fall within the domain of validity of the NRSur7dq4 model, which we also confirm after computing posteriors with the NRSur7dq4 model. A total of 47 BBH events that satisfy these criteria will be used in our analysis; these are listed in Figure 1.
We summarize our parameter inference results using NRSur7dq4 and NRSur7dq4Remnant for the 47 events selected from GWTC-3 [5,6,76,77] according to the procedure laid out in Sec. III. We then compare our results against the public LVK posteriors [5,6,76,77]. While GWTC-2 [4] also provides analyses of a subset of these events, including some results obtained using NRSur7dq4, to compare against the largest possible sample with the most consistent settings, we elide detailed comparisons to this earlier work for simplicity. For similar reasons, we also omit detailed comparison against other analyses of these events [8][9][10][11].
For each event we infer all 15 BBH source parameters discussed in Sec II A, we present resultsfoot_0 for a smaller set of summary observables: the source-frame total mass, mass ratio, dimensionless spin magnitudes χ i , dimensionless spin tilts θ i between the spin vectors and the orbital angular momentum, the effective inspiral spin parameter [103][104][105],
and the transverse spin precession parameter [106,107],
In Fig. 1, we show the recovery of source-frame masses as well as the distance and inclination angles; in the accompanying Fig. 2, we show the corresponding results for
100 200 M [M ] GW150914 GW170729 GW170809 GW170818 GW170823 GW190413 GW190413 GW190421 GW190426 GW190503 GW190513 GW190514 GW190517 GW190519 GW190521 GW190521 GW190527 GW190602 GW190620 GW190630 GW190701 GW190706 GW190727 GW190731 GW190803 GW190805 GW190828 GW190910 GW190915 GW190916 GW190926 GW190929 GW191109 GW191222 GW191230 GW200112 GW200128 GW200129 GW200208 GW200209 GW200216 GW200219 GW200220 GW200220 GW200224 GW200302 GW200311 0.00 0.25 0.50 0.75 1.00 q 50 100 150 200 m 1 [M ]
50 100 m 2 [M ] 10 2 10 3 10 4 D L [Mpc] -1.0 -0.5 0.0 0.5 1.0 cosθ JN
Figure 1. Posteriors for the source-frame total mass M , mass ratio q, source-frame component masses m1, m2, luminosity distance D L and the cosine of the inclination angle θ JN for all 47 events analyzed with NRSur7dq4 model (blue). For comparison, we also show the public LVK posteriors [5,6,76,77] obtained using IMRPhenomXPHM (orange) and SEOBNRv4PHM (green) models. For some events, the SEOBNRv4PHM posteriors are missing from the LVK release (see Tab. I); the SEOBNRv4PHM results are, therefore, also absent in this and all following figures for that subset of events. The grey dashed line represents the mass ratio cut of q = 1/6 used for NRSur7dq4. Posteriors are reported in the wave frame (see Sec. II C) at f ref = 20 Hz. Further details are given in Sec. IV.
0.00 0.25 0.50 0.75 1.00
0.00 0.25 0.50 0.75 1.00 χ 1 0.00 0.25 0.50 0.75 1.00 χ 2
IMRPhenomXPHM SEOBNRv4PHM
Figure 2. Posteriors for the spin magnitudes χ1 and χ2, spin angles θ1 and θ2, effective inspiral spin parameter χ eff , and spin precession parameter χp for all 47 events analyzed with NRSur7dq4 model (blue). For comparison, we also show the public LVK posteriors [5,6,76,77] obtained using IMRPhenomXPHM in orange and SEOBNRv4PHM (where available) in green. Posteriors are reported in the wave frame (see Sec. II C) at f ref = 20 Hz. We provide 3D visualizations of the full spin posteriors at Ref. [78]. Further details are given in Sec. IV.
q [JSD (bits)] IMRPhenomXPHM vs NRSur7dq4 SEOBNRv4PHM vs NRSur7dq4
Figure 3. Jensen-Shannon divergence (JSD) values between the one-dimensional marginalized posteriors of the source-frame total mass M , mass ratio q, source-frame component masses m1, m2, luminosity distance D L and the inclination angle θ JN obtained using NRSur7dq4 and the public LVK posteriors [5,6,76,77] obtained using IMRPhenomXPHM (blue circles) and SEOBNRv4PHM (green squares, where available). Dashed red lines correspond to a JS divergence of 0.02, indicating significant differences between these posteriors. Further details are discussed in Sec. IV. IMRPhenomXPHM vs NRSur7dq4 SEOBNRv4PHM vs NRSur7dq4
Figure 4. Jensen-Shannon divergence (JSD) values between the one-dimensional marginalized posteriors of the spin magnitudes χ1 and χ2, spin angles θ1 and θ2, effective inspiral spin parameter χ eff and spin precession parameter χp obtained using NRSur7dq4 and the public LVK posterior samples [5,6,76,77] obtained using IMRPhenomXPHM (blue circles) and SEOBNRv4PHM (green squares, where available). Dashed red lines correspond to a JS divergence of 0.02, indicating significant differences between these posteriors. Further details are discussed in Sec. IV.
the inferred spin parameters. For comparison, in both figures we also show the results from the public LVK posteriors [5,6,76,77] obtained using the IMRPhenomXPHM and SEOBNRv4PHM modelsfoot_1 . Note that IMRPhenomXPHM posteriors are obtained with bilby [99] while SEOBNRv4PHM posteriors are obtained with RIFT [108]. For the events listed in Table I, SEOBNRv4PHM results are absent from the public LVK data release [5,6,76,77] and, therefore, are also absent in our comparisons. In Figs. 1 and 2, we find many events for which noticeable differences exist between posteriors inferred with NRSur7dq4, IMRPhenomXPHM, and SEOBNRv4PHM. We use a standard diagnostic -the Jensen-Shannon (JS) divergence [109] -to quantify the difference between the one-dimensional marginalized posteriors inferred with NRSur7dq4 model (this work) and the posteriors obtained using IMRPhenomXPHM and SEOBNRv4PHM models (from the publicly available LVK posteriors [5,6,76,77]). Recall the JS divergence (JSD) between two probability density functions p(x) and q(x),
is defined as a symmetrized extension of the Kullback-Leibler divergence [110], where m(x) = 1/2(p + q) is the point-wise mean of p(x) and q(x) and KLD,
is Kullback-Leibler divergence. When using a base-2 logarithm JS divergence values are given in the units of bits. A JS divergence value of 0 bits signifies that the posteriors are identical, while a JS divergence value of 1 bit corresponds to posterior distributions that have no statistical overlap at all. Values smaller than 2 × 10 -3 bits can occur due to stochastic sampling and indicate statistically indistinguishable posterior samples [79]. The threshold for non-negligible bias varies by study, where JSD values of 0.007 [4] (which for a Gaussian corresponds to a 20% shift in the mean), 0.02 [79], 0.05 [111], and 0.15 [112] have all been used. In this paper, we will consider JSD values above 0.02 bits to indicate important differences between posteriors recovered from two different BBH waveform models [79]. In Figs. 34, we show the JS divergence values between posterior distributions for the masses, spin, distance, and inclination angles (i.e., for the parameters shown in Figs. 12), for each event. We find that the JS divergence values between NRSur7dq4 and the IMRPhenomXPHM and SEOBNRv4PHM models are mostly less than 0.02 bits suggesting good agreement. However, for around ∼ 23% (∼ 55%) of the analyzed events, JS divergence values between NRSur7dq4 and IMRPhenomXPHM (SEOBNRv4PHM) are larger than 0.02 bits for at least one of the parameters shown in Figs. 34. Such differences can arise from waveform systematics, such as the missing physics in IMRPhenomXPHM/SEOBNRv4PHM from not being informed by precessing NR simulations. Some of the observed differences may stem from the sampler techniques used:
bilby was employed for IMRPhenomXPHM, RIFT for SEOBNRv4PHM, and parallel-bilby for NRSur7dq4. Differences between samplers are expected to be most significant when comparing results obtained with RIFT to those from parallel-bilby, as their underlying methods differ substantially. In contrast, parallel-bilby and bilby share the same nested sampling algorithm, leading to greater consistency between their results. In Appendix A, we validate the agreement between our parallel-bilby parameter estimation setup and the LVK's bilby setup. Therefore, any discrepancies in inference results between IMRPhenomXPHM and NRSur7dq4 are likely due to modeling systematics.
To ensure that sufficient sampling density exists in the publicly available SEOBNRv4PHM posteriors computed by RIFT, we confirmed that the Jensen-Shannon divergence (JSD) values reported in Figs. 3 and 4 remain reliable even for events with a limited number of posterior samples. As a test, for each event we downsampled the SEOBNRv4PHM posterior by 80% and recalculated all JSD values. The resulting JSD values changed by less than 0.01, with most differences being below 0.001. We further checked that the recomputed JSD values are visually indistinguishable from the results already presented in the figures. Thus, we conclude that the reported JSD values are robust and reliable, even for the apparently undersampled SEOBNRv4PHM posteriors. Further investigation, however, is needed to completely identify the specific sources of differences between posteriors computed with SEOBNRv4PHM and NRSur7dq4.
Having provided an overview of the comparison of NRSur7dq4 posteriors against SEOBNRv4PHM and IMRPhenomXPHM results in Sec. IV, we now highlight certain events with noticeable differences in their inferred posteriors. For these events, Fig. 5 shows the posteriors for mass and spin parameters, while Fig. 6 focuses on the distance and inclination posteriors (with the total mass shown for comparison). Median values of the inferred source properties are given in Table III in the Appendix.
We reiterate that, for some of the events, SEOBNRv4PHM posteriors (obtained using the RIFT parameter estimation code) are missing in public LVK posteriors, as listed in Table I. Furthermore, for several events where SEOBNRv4PHM posteriors are available due to an inefficient post-processing step incorporating the effects of calibration uncertainties in the GW data, the number of samples can be significantly reduced [113,114]. This results in the under-sampled posteriors for SEOBNRv4PHM visible for some events shown in Fig. 5 and Fig. 6.
GW150914_095045 was the first direct GW observation and, with an SNR of ∼ 26 [111], remains one of the loudest events detected. This event has been extensively studied using a variety of waveform models over the past several years [70,111,115,116]. We find that while most of the marginalized posterior distributions match across waveform models, the distributions describing the spin parameters do not, as can be seen in the top row of Fig. 5.
In Fig. 5 (top row), we show that NRSur7dq4 favors smaller values for the spin magnitudes χ 1 and χ 2 compared to IMRPhenomXPHM and SEOBNRv4PHM. We also observe significant differences in the χ p posteriors (Fig. 7a), where NRSur7dq4 favors smaller χ p values than IMRPhenomXPHM and SEOBNRv4PHM. Fig. 7b summarizes the JS divergence values for many of the most interesting parameters. Interestingly, as highlighted in Fig. 7a, the NRSur7dq4 posterior for χ p matches more closely with earlier estimates from GWTC-1 [111] using the IMRPhenomPv2 [38, 39, 106] (which does not include higher order modes) and SEOBNRv3 [117][118][119] (which do not include ℓ > 2 modes) models. While we cannot offer a simple explanation for this tension between models, especially as the different models also rely on different posterior samplers 5 , this might point to systematic differences arising from the treatment of precession or whenever subdominant harmonic modes play a role.
For GW190413_13430 we find noticeable differences in the mass ratio (second row of Fig. 5), luminosity distance and inclination (first column of Fig. 6). The JS divergence values between the one-dimensional marginalized posteriors recovered with NRSur7dq4 and IMRPhenomXPHM are 0.031 bits (for q), 0.05 bits (for D L ), and 0.026 bits (for θ JN ). We cannot compare with SEOBNRv4PHM because the corresponding posteriors are missing from the most recent LVK release.
GW190521_030229 [101,121] has unusually high component masses, 85 +21 -14 M ⊙ and 66 +17 -18 M ⊙ , compared to other events. Consequently, the observable signal contains only a few pre-merger cycles. Furthermore, this event shows tentative signs of precession [101,[121][122][123][124]. As a result, systematic differences between waveform models in the treatment of both precession and the merger-ringdown can play a more significant role for signals like this one. For example, Refs. [73,123,124] already pointed out several systematic differences between NRSur7dq4 and other models for this event, and the same can be seen in Figs. 5 (third row) and 6 (second column). In particular, IMRPhenomXPHM prefers more unequal mass ratios with a clear bimodal mass ratio distribution while the NRSur7dq4 posterior shows visibly different features [123,124] (third row of Fig. 5). Similarly, the posteriors for several spin parameters, including χ 1 , χ 2 , θ 1 , θ 2 , χ eff and χ p in Fig. 5 indicate significant systematic differences between NRSur7dq4 and IMRPhenomXPHM. The JS divergence values between the one-dimensional marginalized posteriors recovered with NRSur7dq4 and IMRPhenomXPHM are 0.043 bits (for M ), 0.136 bits (for q), 0.002 bits (for χ 1 ), 0.025 bits (for χ 2 ), 0.102 bits (for θ 1 ), 0.017 bits (for θ 2 ), 0.126 bits (for χ eff ) and 0.092 bits (for χ p ). Among the extrinsic parameters, we find significant differences between the posteriors for D L and θ JN (second column of Fig. 6), with JS divergence values of 0.083 bits and 0.107 bits, respectively. While we cannot compare with SEOBNRv4PHM as the corresponding posteriors are missing from the most recent LVK release [5,6,76,77], such a comparison is shown in Ref. [73], where significant differences were also noted between SEOBNRv4PHM and NRSur7dq4.
Figures 5 (fourth row) and 6 (third column) show noticeable differences between NRSur7dq4 and the other models for several parameters for this event. In particular, NRSur7dq4 and SEOBNRv4PHM yield similar posteriors for the source-frame total mass M whereas IMRPhenomXPHM favors slightly larger values for M . NRSur7dq4 also favors a slightly more asymmetric binary i.e., smaller values for q compared to the other two models. Furthermore, we find noticeable differences in the χ 1 , χ 2 , and χ p posteriors between three models, with the NRSur7dq4 falling broadly in between IMRPhenomXPHM and SEOBNRv4PHM. The JS divergence values between NRSur7dq4 and IMRPhenomXPHM (SEOBNRv4PHM) posteriors for M , q, χ 1 , χ 2 and χ p are 0.04 bits, 0.016 bits, 0.036 bits, 0.013 bits, and 0.047 bits (0.004 bits, 0.035 bits, 0.042 bits, 0.005 bits, and 0.036 bits), respectively. We further find noticeable differences between NRSur7dq4 and IMRPhenomXPHM (SEOBNRv4PHM [to a lesser extent]) posteriors for D L and θ JN (third column of Fig. 6), with JS divergence values of 0.172 bits and 0.234 bits (0.04 bits and 0.07 bits), respectively.
58.0 64.5 71.0 M [M ] GW150914 095045 0.1 0.55 1.0 q 0.0 0.5 1.0 χ1 0.0 0.5 1.0 χ2 0.0 π/2 π θ1 [rad] 0.0 π/2 π θ2 [rad] -0.5 0.0 0.5 χeff 0.0 0.5 1.0 χp 60.0 91.0 122.0 M [M ] GW190413 134308 0.1 0.55 1.0 q 0.0 0.5 1.0 χ1 0.0 0.5 1.0 χ2 0.0 π/2 π θ1 [rad] 0.0 π/2 π θ2 [rad] -0.7 0.0 0.7 χeff 0.0 0.5 1.0 χp 122.0 199.0 276.0 M [M ] GW190521 030229 0.1 0.55 1.0 q 0.0 0.5 1.0 χ1 0.0 0.5 1.0 χ2 0.0 π/2 π θ1 [rad] 0.0 π/2 π θ2 [rad] -0.8 0.0 0.8 χeff 0.0 0.5 1.0 χp 67.0 78.5 90.0 M [M ] GW190521 074359 0.1 0.55 1.0 q 0.0 0.5 1.0 χ1 0.0 0.5 1.0 χ2 0.0 π/2 π θ1 [rad] 0.0 π/2 π θ2 [rad] -0.4 0.0 0.4 χeff 0.0 0.5 1.0 χp 42.0 97.5 153.0 M [M ] GW190527 092055 0.1 0.55 1.0 q 0.0 0.5 1.0 χ1 0.0 0.5 1.0 χ2 0.0 π/2 π θ1 [rad] 0.0 π/2 π θ2 [rad] -0.7 0.0 0.7 χeff 0.0 0.5 1.0 χp 85.0 122.0 159.0 M [M ] GW191109 010717 0.1 0.55 1.0 q 0.0 0.5 1.0 χ1 0.0 0.5 1.0 χ2 0.0 π/2 π θ1 [rad] 0.0 π/2 π θ2 [rad] -0.8 0.0 0.8 χeff 0.0 0.5 1.0 χp 57.0 65.0 73.0 M [M ] GW200129 065458 0.1 0.55 1.0 q 0.0 0.5 1.0 χ1 0.0 0.5 1.0 χ2 0.0 π/2 π θ1 [rad] 0.0 π/2 π θ2 [rad] -0.3 0.0 0.3 χeff 0.0 0.5 1.0 χp NRSur7dq4 IMRPhenomXPHM SEOBNRv4PHM
Figure 5. Posteriors for the source-frame total mass M , mass ratio q, spin magnitudes χ1 and χ2, spin angles θ1 and θ2, effective inspiral spin parameter χ eff , and spin precession parameter χp for a list of seven events for which we infer the most significant differences between results obtained using NRSur7dq4 (blue histogram), IMRPhenomXPHM (orange histogram) and SEOBNRv4PHM (green histogram, where available). The under-sampled posteriors for SEOBNRv4PHM is a consequence of an inefficient post-processing procedure used in the RIFT code whenever calibration uncertainties are accounted for [113,114]. Posteriors are reported in the wave frame (see Sec. II C) at f ref = 20 Hz. Further details are given in Sec. IV A.
60.0 91.0 122.0 M [M ] GW190413 134308 1000.0 5000.0 9000.0 DL [Mpc] -1.0 0.0 1.0 cosθJN 122.0 199.0 276.0 M [M ] GW190521 030229 1000.0 5000.0 9000.0 DL [Mpc] -1.0 0.0 1.0 cosθJN 67.0 78.5 90.0 M [M ] GW190521 074359 300.0 1150.0 2000.0 DL [Mpc] -1.0 0.0 1.0 cosθJN 85.0 122.0 159.0 M [M ] GW191109 010717 350.0 1925.0 3500.0 DL [Mpc] -1.0 0.0 1.0 cosθJN 57.0 65.0 73.0 M [M ] GW200129 065458 350.0 875.0 1400.0 DL [Mpc] -1.0 0.0 1.0 cosθJN NRSur7dq4 IMRPhenomXPHM SEOBNRv4PHM
Figure 6. Posteriors for the source-frame total mass M , luminosity distance D L , and cosine of the inclination angle θ JN for a list of five events for which we see the largest differences between results obtained using NRSur7dq4 (blue histogram), IMRPhenomXPHM (orange histogram) and SEOBNRv4PHM (green histogram, where available). Further details are given in Sec. IV A.
GW190527_092055 presents an interesting case as we find that, for almost all parameters shown in Fig. 5 (fifth row), NRSur7dq4 and IMRPhenomXPHM yield consistent posteriors while SEOBNRv4PHM posteriors show noticeable differences. For example, both NRSur7dq4 and IMRPhenomXPHM posteriors for spin magnitudes χ 1 and χ 2 are uninformative whereas SEOBNRv4PHM posteriors show strong support for smaller values of χ 1 and χ 2 (fifth row of Fig. 5). The JS divergence values between NRSur7dq4 and IMRPhenomXPHM (SEOBNRv4PHM) posteriors for χ 1 and χ 2 are 0.0004 bits and 0.0003 bits (0.05 bits and 0.03 bits) respectively. Posteriors for the spin angles θ 1 and θ 2 , however, match for all models (fifth row of Fig. 5). However, the SEOBNRv4PHM posteriors (obtained using the RIFT code) for this event appear to be particularly undersampled, making it difficult to disentangle model systematics from sampler systematics. We note that Ref. [125] reanalyzed this event using parallel-bilby with NRSur7dq4 and obtained results consistent with ours. However, Ref. [126] employed a machine-learning based parameter estimation code [127] with importance sampling to reanalyze this event with SEOBNRv4PHM, finding better agreement between SEOBNRv4PHM and IMRPhenomXPHM.
GW191109_010717 is another event that shows interesting and astrophysically important differences between posteriors obtained using the NRSur7dq4, IMRPhenomXPHM, and SEOBNRv4PHM models. For mass ratio q and spin magnitude χ 1 , IMRPhenomXPHM posteriors show bimodalities whereas NRSur7dq4 and SEOBNRv4PHM posteriors do not (sixth row of Fig. 5). Another noteworthy observation is that IMRPhenomXPHM favors larger values for the secondary spin magnitude χ 2 (sixth row of Fig. 5), whereas both NRSur7dq4 and SEOBNRv4PHM present posteriors for χ 2 that are effectively uninformative. Furthermore, NRSur7dq4 shows a stronger preference for negative χ eff at 99.3% credible level, compared to 95.9% for SEOBNRv4PHM 6 and 85.3% for IMRPhenomXPHM. For IMRPhenomXPHM, we also see a bimodality in χ eff , which likely results from the bimodality in q, due to the correlation between these two parameters. Furthermore, we note that the spin angle θ 1 is well measured for this event but noticeably different across waveform models. Finally, the χ p posteriors are also noticeably different across the models. The JS divergence values between NRSur7dq4 and IMRPhenomXPHM (SEOBNRv4PHM) posteriors for q, χ 1 , χ 2 , χ eff , χ p and θ 1 are 0.091 bits, 0.062 bits, 0.067 bits, 0.139 6 Interestingly, Ref. [62] found that the newer SEOBNRv5PHM model shows a stronger preference for negative χ eff for GW191109_010717, in agreement with NRSur7dq4.
bits, 0.029 bits and 0.082 bits (0.072 bits, 0.012 bits, 0.009 bits, 0.086 bits, 0.056 bits and 0.117 bits) respectively. Among the extrinsic parameters, NRSur7dq4 provides more tightly constrained posteriors for the luminosity distance D L and inclination θ JN (fourth column of Fig. 6). The JS divergence values between NRSur7dq4 and IMRPhenomXPHM (SEOBNRv4PHM) posteriors for D L and θ JN are 0.081 bits and 0.107 bits (0.278 bits and 0.138 bits) respectively.
The strong preference for χ eff < 0 can have important astrophysical implications, as negative χ eff is expected to be more common in dynamically formed binaries than those formed through isolated evolution [128,129]. We note, however, that de-glitched strain data is used for this event. Previous work on GW200129_065458 has shown potentially subtle issues can arise when using deglitched strain data [130]. We refer to the dedicated glitch subtraction study presented for this and several other events in Ref. [131] and Ref. [132]. Importantly, Ref. [132] found (see their App. A) that transient nongaussian noise or glitches affecting the data around the time of this event led to false violations of GR in the tests conducted in that work. Such effects could also impact the inference of χ eff (see e.g. App. B of Ref. [131]).
Next, we look at GW200129_065458 -another event where glitch subtraction was necessary [90,130], thereby complicating a straightforward interpretation of the signal. This event has many interesting properties: it has a network matched-filter SNR of 26.8 [6] making it the loudest detected BBH signal, it has observable spin-induced orbital precession [75], and the post-merger remnant BH has a large recoil velocity [74].
Comparing the posteriors for NRSur7dq4, SEOBNRv4PHM, and IMRPhenomXPHM, we find noticeable differences for several parameters in Fig. 5 (seventh row) and Fig. 6 (fifth column). For example, NRSur7dq4 and IMRPhenomXPHM posteriors exhibit varying degrees of bimodality in mass ratio q while SEOBNRv4PHM posteriors are uni-modal. Furthermore, NRSur7dq4 and IMRPhenomXPHM favor smaller (more unequal) values of q than SEOBNRv4PHM. We also find significant differences in χ p estimates between NRSur7dq4 and SEOBNRv4PHM (seventh row of Fig. 5) while IMRPhenomXPHM posteriors are consistent with NRSur7dq4. We note that our NRSur7dq4 posteriors for χ p match the results obtained in Ref. [75]. Similarly, for θ 1 and θ 2 , NRSur7dq4 and SEOBNRv4PHM posteriors show noticeable differences whereas IMRPhenomXPHM broadly agrees with NRSur7dq4. The JS divergence values between NRSur7dq4 and IMRPhenomXPHM (SEOBNRv4PHM) posteriors for q, χ 1 , χ 2 , χ eff , χ p , θ 1 and θ 2 are 0.023 bits, 0.002 bits, 0.017 bits, 0.063 bits, 0.001 bits, 0.008 bits and 0.008 bits (0.378 bits, 0.331 bits, 0.033 bits, 0.129 bits, 0.405 bits, 0.174 bits and 0.09 bits) respectively.
Finally, we also find significant differences between For comparison, we also show posteriors obtained in earlier LVK studies (GWTC-1) using IMRPhenomPv2 (light orange) and SEOBNRv3 (light green). Somewhat surprisingly, the NRSur7dq4 posterior agrees more closely with older IMRPhenomPv2 and SEOBNRv3 models that do not include ℓ > 2 subdominant modes. (b) Jensen-Shannon divergence (JSD) values between the one-dimensional marginalized posteriors for a set of parameters (shown in Fig. 3 and Fig. 4) obtained using NRSur7dq4 and the public LVK posterior samples [5,6,76,77] obtained using IMRPhenomXPHM (blue circles) and SEOBNRv4PHM (green squares). Dashed red lines correspond to a JS divergence of 0.02, indicating significant differences between these posteriors. Further details are given in Sec. IV A 1.
NRSur7dq4 and IMRPhenomXPHM (SEOBNRv4PHM) posteriors for D L and θ JN (fifth column of Fig. 6), with JS divergence values of 0.025 bits and 0.029 bits (0.124 bits and 0.228 bits), respectively. We note that NRSur7dq4 favors larger values for D L than both the IMRPhenomXPHM and SEOBNRv4PHM models.
One generally expects the JS divergence between posteriors inferred using different waveform models to be larger for higher-SNR signals and smaller for lower-SNR signals. This follows from the fact that lower-SNR events provide less informative data, leading to posteriors that remain closer to the prior distribution. Additionally, certain parameters -such as the secondary spin, χ 2 -are known to be more difficult to measure, resulting in typically lower JS divergence values compared to parameters like the total mass or effective spin, which tend to be better constrained.
We examine the dependence of the JS divergence on the network SNR for posteriors obtained using the NRSur7dq4 model. Figure 8a confirms the expected trend: higher-SNR events generally exhibit greater differences between posteriors inferred with the NRSur7dq4 and IMRPhenomXPHM models. Notably, five of the seven highlighted events from Sec. IV A have SNRs exceeding 14, placing them in the upper tail of the SNR distribution. However, some high-SNR events deviate from this pattern, suggesting that in certain regions of the binary black hole parameter space, the NRSur7dq4 and IMRPhenomXPHM models exhibit relatively low systematic differences.
A similar trend is observed when comparing NRSur7dq4 and SEOBNRv4PHM, as shown in Fig. 8b. One particularly notable outlier is GW190527_092055, which has a relatively low SNR of about 8 yet exhibits a large JS divergence between the NRSur7dq4 and SEOBNRv4PHM posteriors. As discussed in Sec. IV A, the posteriors for this event, obtained using the RIFT code, may be dominated by sampler effects.
Next, we investigate whether NRSur7dq4 posteriors yield better constrained sky localization when com- 1 and 2) between posteriors inferred using the NRSur7dq4 and IMRPhenomXPHM models as a function of network matched-filter SNR. As expected, higher-SNR events generally exhibit larger JS divergence, indicating greater sensitivity to waveform model differences. Events highlighted in orange correspond to those discussed in Sec. IV A. While most high-SNR events follow the expected trend, some deviate, suggesting regions of the binary black hole parameter space where these two models show less systematic disagreement. (b) Same as Fig. 8a, but comparing posteriors obtained with the NRSur7dq4 and SEOBNRv4PHM models. The general trend remains, with larger JS divergence for higher-SNR events. Notably, however, GW190527_092055 (SNR ≈ 8) stands out as an outlier, exhibiting a relatively large JS divergence. This event's posteriors, computed using the RIFT code, may be dominated by sampler effects.
pared to the public LVK posteriors obtained using the IMRPhenomXPHM and SEOBNRv4PHM models [5,6,76,77]. For most events, NRSur7dq4, IMRPhenomXPHM, and SEOBNRv4PHM offer largely consistent skymaps.
In Fig. 9, we show the recovered skymaps for three events for which we notice the largest differences; the contours show regions containing the central 50% and 90% of the two-dimensional posterior distribution over sky angles -right ascension α and declination δ. For each of the three events considered in Fig. 9, we find that the skymaps for NRSur7dq4 and IMRPhenomXPHM are consistent, while SEOBNRv4PHM shows a significant difference. It is important to recall that NRSur7dq4 and IMRPhenomXPHM posteriors are computed using parallel-bilby and bilby, respectively, and both employ the same dynesty sampler. The SEOBNRv4PHM posteriors are obtained with RIFT, which employs different sampling techniques over the extrinsic parameters. This is a potential reason for the differences in the skymap posteriors of SEOBNRv4PHM compared to the other models.
To understand whether the data prefers a particular waveform model, one can compare the Bayes factors, given in Eq. ( 6), for different models. For simplicity, we assume all models are equally likely, thereby sidestepping the issue of setting prior model odds. Even with this simplification, meaningfully comparing Bayes factors is complicated by the fact that the prior used for NRSur7dq4 in our study is restricted to a smaller portion of the parameter spacefoot_3 as compared to the ones used for IMRPhenomXPHM and SEOBNRv4PHM in the LVK analyses [5,6,76,77]. Therefore, we re-analyze all 47 events considered in this paper using the IMRPhenomXPHM model with the same restricted priors and sampler settings (see Sec.II F) used for the NRSur7dq4 runs. We do not perform any new parameter estimation runs with SEOBNRv4PHM due to the model's high computational cost. Note that these IMRPhenomXPHM results obtained with redistricted priors are used only for
Figure 9. Skymaps for three events where noticeable differences are observed between NRSur7dq4 (blue), IMRPhenomXPHM (orange), and SEOBNRv4PHM (green). We note that the discrepancies between SEOBNRv4PHM and the other models may arise from a different sampler used for the analysis with SEOBNRv4PHM.
Further details are given in Sec IV C.
meaningfully comparing the Bayes factors and SNRs in Fig. 10 and Fig. 18. For all other comparisons made throughout this paper, we use IMRPhenomXPHM results from the public LVK posteriors [5,6,76,77].
We compute the differences,
∆ log e B = log e B NRSurlog e B XPHM , (11) where Eq. ( 11) is motivated by the identity
Here B NRSur is the Bayes factor of NRSur7dq4 over noise hypothesis and B XPHM is the Bayes factor of IMRPhenomXPHM over noise hypothesis using the same priors used for the NRSur7dq4 analyses; see Eq. ( 6).
The log Bayes factor differences between NRSur7dq4 and IMRPhenomXPHM is shown in the upper panel of Fig. 10. We further provide associated error estimates for ∆ log e B (shaded red region), obtained using the error estimates for the logarithm of the waveform-model evidence provided by dynesty, and adding them in quadrature. We note, however, that this error band provides only a rough guide to indicate the accuracy achieved when computing the Bayes factor [97,[133][134][135]. Indeed, we should view dynesty's method for computing Bayes factors probabilistically, and the error estimation as a probabilistic statement (say, 1-sigma interval) instead of an error bound. For instance, in Fig. 10, we compare dynesty's error estimate with the difference between the maximum and minimum log-Bayes factors (shaded green region) from the four independent Bayesian inference runs with different random seeds performed for each event (see Sec. II F). This provides an an alternative estimate of the errors in the log-Bayes factors 8 . We find that, for at least some events, the difference in log Bayes factor between the runs with different seeds can be larger than the error estimate in log Bayes factor provided by dynesty/parallel-bilby. This is not surprising in light of dynesty's probabilistic error estimation [97,[133][134][135], but it does mean that for any particular event a log Bayes factor value can fluctuate outside of the red shaded region (dynesty error estimate) due to stochastic sampling. The histogram's bin size has been set to roughly a 1-sigma interval.
Combining all 47 events in Fig. 10, the cumulative ∆ log e B value for NRSur7dq4 over IMRPhenomXPHM is 0.54, suggesting an overall mild preference for NRSur7dq4. We find three events, GW191109_010717, GW190521_030229 and GW190521_074359, where there is a clear preference for NRSur7dq4, with ∆ log e B values of 2.73, 1.69, and 1.08, respectively. All three of these events were highlighted in Sec. IV A as showing noticeable differences in posteriors. Excluding these three events, the histogram in the top-right panel of Fig. 10 indicates that (i) most events show no preference for either model, and (ii) many events show a very mild preference for IMRPhenomXPHM, although, as noted above, the true errors in our Bayes factor computation may be larger than those indicated in Fig. 10 thereby spoiling a clear interpretation of small values for any particular between NRSur7dq4 (referred to as 'NRSur') and IMRPhenomXPHM models (referred to as 'XPHM') for all 47 events. These events were analyzed using the prior described in Sec. II E, which is restricted to a smaller portion of the parameter space as compared to the official LVK prior but matches the one used for our NRSur7dq4 analysis. We also show the estimated uncertainty in the calculation of the log evidence from parallel-bilby as a shaded red patch (upper panel). On the other hand, the shaded green patch (upper panel) shows the uncertainties estimated from four independent Bayesian inference runs with different random seeds. Further details are given in Sec. V.
event. GW190706_222641 shows the largest preference for IMRPhenomXPHM, with a ∆ log e B value of -0.98.
In summary, while the data shows a very mild preference for NRSur7dq4 over IMRPhenomXPHM, there are clear outlier events and secular trends in the distribution of Bayes factors. Outlier events, trends, and the robustness of our results are considered in App. B.
Next, we compute posteriors for the network-matched filter SNR, ρ, recovered by the NRSur7dq4 model for each event. We then compare the median values of the posteriors for ρ against the ones obtained using IMRPhenomXPHM. In the lower panel of Fig. 10, we report the difference,
between the median network matched filter SNR recovered by NRSur7dq4 and IMRPhenomXPHM. In this case, as indicated in the histogram in the bottom-right panel of Fig. 10, we find NRSur7dq4 typically recovers larger SNRs. GW191109_010717 and GW190521_030229 show the largest preference for NRSur7dq4, with differences in median SNR of 0.28 and 0.22, respectively. On the other hand, GW190517_055101 shows the largest preference for IMRPhenomXPHM, with a difference in median SNR of -0.20. Outlier events, trends, and the robustness of our results are considered in App. B.
In addition to constraining the properties of the component BHs in the binary, we infer the properties of the final BH left behind after the merger, in particular, its source-frame mass M f , spin vector χ f , and recoil velocity vector v f . During its evolution, the binary radiates energy, angular momentum, and linear momentum. The radiated energy and angular momentum are reflected in M f and χ f , respectively. The radiated linear momentum
50 150 250 0.25 0.50 0.75 1.00
. Posteriors for the source-frame remnant mass M f , the remnant spin magnitude χ f , the remnant kick magnitude v f , and the spherical polar {θχ f , θv f } and azimuthal angles {ϕχ f , ϕv f } of χ f and v f for all 47 events analyzed with the NRSur7dq4 model (blue). For comparison, we also show the M f and χ f posteriors from the public LVK data release [5,6,76,77] obtained using IMRPhenomXPHM (orange) and SEOBNRv4PHM (green, where available). For the rest of the parameters, we show the effective priors (in gray); the difference between the prior and posterior can be used to assess how informative the data are about these parameters. The spin and kick angles are shown in the wave frame at t ref = -100 M det . We provide 3D visualizations of the full remnant spin and kick posteriors at Ref. [78]. Further details are given in Sec. VI.
causes a shift in the binary's center of mass in the opposite direction, imparting a recoil velocity or a "kick" v f to the remnant BH. While χ f is restricted to be either parallel or anti-parallel to the orbital angular momentum L for nonprecessing binaries, it can be arbitrarily oriented for precessing binaries. On the other hand, while the kicks for nonprecessing binaries are typically restricted to ≲ 300 km/s and along the orbital plane, kicks for precessing binaries can reach magnitudes up to ∼ 5000 km/s [136][137][138][139] with arbitrary orientations. However, as we will discuss below, for precessing binaries, the direction of χ f is preferentially along, while the direction of v f is preferentially along or opposite, the direction of L near merger.
Remnant BH properties have important applications for astrophysics and fundamental physics. The remnant mass and spin magnitude are important for tests of general relativity using GWs [132], as the the remnant mass entirely determines frequencies in the ringdown and spin. The kick magnitude is important for placing observational constraints [19,73,74,[140][141][142] on the rate of hierarchical mergers in dense environments: repeated mergers are a means to form heavy BHs in nature, but if the kick exceeds the escape velocity of the host environment, the remnant BH after the first merger would simply get ejected and not participate in another merger. Finally, the remnant spin direction and kick direction [74,143] can be useful to study binaries that may be formed in active galactic nuclei disks, as the final BH's spin orientation with respect to the disk as well as its motion can impact potential electromagnetic counterparts [144]. The kick direction also shows up when computing the Dopplershifted remnant mass, which may play a role in future high-accuracy ringdown tests of GR [74,140,145,146].
Therefore, while the public LVK results [5,6,76,77] only include the remnant mass and spin magnitudes, we provide posterior samples for the full spin and kick vectors. We report the source-frame mass M f , spin vector χ f and kick velocity vector v f of the remnant black hole for each of the 47 events we analyze in Fig. 11.
We follow the prescription outlined in Refs. [67,74,86,140]: starting with the posterior samples for the spins and detector-frame component masses obtained using the NRSur7dq4 model, we evaluate the associated remnant surrogate model NRSur7dq4Remnant which provides estimates for M f,det , χ f and v f , from which we compute M f = M f,det /(1 + z). These samples are included in our public release at Ref. [78]. Note that the remnant spin and kick vectors in our public release are defined in the same frame as the component BHs, i.e. the wave frame at f ref = 20 Hz, as described in Sec. II C.
However, when visualizing the remnant spin and kick directions in this section, we adopt a different frame that is more naturally suited for discussing remnant properties: the wave frame at t ref = -100 M det as proposed in Ref. [72]. This frame is similar to the wave frame at f ref = 20 Hz (see Sec. II C), except that the reference point is chosen to be the dimensionless time of t ref = -100 M det before the peak waveform amplitude (defined in Eq.5 of Ref. [67]). Because this reference point is always very close to the merger (typically within 2-4 GW cycles [72]), it provides a more natural frame to define remnant spin and kick vectors than the wave frame at f ref = 20 Hz (which can occur up to ∼ 40 GW cycles before the mergers for the events considered in this work).
For example, in the wave frame at t ref = -100 M det , the direction of χ f is preferentially oriented close to the z-axis. This can be explained as follows: the direction of χ f can be approximated by the direction [147,148] of the total angular momentum J = L + m 2 1 χ 1 + m 2 2 χ 2 , but J is typically dominated (excluding the special case of transitional precession) by the contribution from L rather than the contribution from the spins. As a result, the direction of χ f is preferentially oriented close to L near merger, which is along the z-axis of the wave frame at t ref = -100 M det (see Sec. II C). This is reflected in the prior for the remnant spin direction, as we will see in Sec. VI D. Similarly, the wave frame at t ref = -100 M det is well-suited to discuss the kick direction as well, as the kick is known to be preferentially orientated close to or opposite to L near merger [86,[136][137][138][139]. For this reason, while our public release [78] will contain remnant spin and kick vectors defined in the wave frame at f ref = 20 Hz for consistency with the frame used for the component BH spins, in Figs. 11, 12, 13, 14, 15, we will adopt the wave frame at t ref = -100 M det . As described in Ref. [72], the two frames are related by a transformation described by the dynamics of NRSur7dq4, which is provided by the model [67].
Figure 11 summarizes our constraints on the remnant properties for all 47 events considered in our analysis (cf. Sec. III). The first two columns show the sourceframe remnant mass M f and spin magnitude χ f for NRSur7dq4 with NRSur7dq4Remnant, along with the corresponding constraints from the LVK public release for IMRPhenomXPHM and SEOBNRv4PHM [5,6,76,77]. In the LVK results, the remnant mass and spin magnitude are computed following Ref. [149], using the remnant models of Refs. [148,150,151]. While these models include some precession corrections, they are not informed by precessing NR simulations. Therefore, differences with respect to our estimates of M f and χ f can arise from differences in the remnant models as well as the waveform model used to infer binary source parameters.
The first two columns of Fig. 12 show the JS divergence between our posteriors for M f and χ f , and the public LVK samples for IMRPhenomXPHM and SEOBNRv4PHM. For about 13 % (17% ) of events, the JS divergence between NRSur7dq4 and IMRPhenomXPHM (SEOBNRv4PHM)
56.0 61.5 67.0 Mf [M ] GW150914 095045 0.5 0.6 0.8 χf 116.0 186.5 257.0 Mf [M ] GW190521 030229 0.2 0.6 0.9 χf 41.0 95.5 150.0 Mf [M ] GW190527 092055 0.2 0.6 0.9 χf 82.0 116.0 150.0 Mf [M ] GW191109 010717 0.2 0.6 0.9 χf 54.0 61.5 69.0 Mf [M ] GW200129 065458 0.6 0.7 0.8 χf NRSur7dq4 IMRPhenomXPHM SEOBNRv4PHM Figure 13. Posteriors of source-frame mass and spin magnitude of the remnant BH for five events where we see the most prominent differences between results obtained using NRSur7dq4 (blue histogram), IMRPhenomXPHM (orange histogram) and SEOBNRv4PHM (green histogram, where available). Further details are given in Sec. VI B.
for M f rises above 0.02 bits, indicating noticeable differences. Similarly, for about 19% (38%) of events, the JS divergence between NRSur7dq4 and IMRPhenomXPHM (SEOBNRv4PHM) for χ f rises above 0.02 bits. The events with the most prominent differences are highlighted in Fig. 13. Constraints on the remaining remnant parameters (the direction of χ f , and the vector v f ) are not provided in the public LVK samples [5,6,76,77]. Therefore, in Figs. 11 and 12, we compare our posteriors for these parameters with their effective priors, to judge how informative the data are about these quantities. Following Refs. [74,140], effective prior samples for M f (not shown), χ f and v f are obtained by evaluating the NRSur7dq4Remnant model on samples drawn from the prior on {m 1 , m 2 , χ 1 , χ 2 } (see Sec. II E). In the following figures, the prior samples for χ f and v f are also transformed to the wave frame at t ref = -100 M det . To parameterize the remnant spin and kick directions, we adopt the standard spherical polar angles θ χ f and θ v f , and azimuthal angles ϕ χ f and ϕ v f , computed in the wave frame at t ref = -100 M det .
The fifth column of Fig. 11 shows our constraints on the kick magnitude v f along with the corresponding prior. As expected from Ref. [140], for most events the posterior is largely indistinguishable from the prior, meaning that the data are not informative about the kick magnitude. In Fig. 12, we find that the JS divergence between the posterior and prior of v f rises above 0.02 bits for about 36% of the events considered. The events with the four highest JS divergence values are highlighted in Fig. 14. Notably, two of these events show a clear preference away from v f = 0 compared to the prior: GW200129_065458 with v f ∼ 1392 +848 -1085 km/s and a JS divergence of 0.305 bits, and GW191109_010717 with v f ∼ 485 +668 -252 km/s and a JS divergence of 0.103 bits. Here, we report the median and 90% symmetric credible interval. While this finding can have important astrophysical implications, especially for hierarchical mergers, we again point out that both GW200129_065458 and GW191109_010717 suffered from being coincident with a detector glitch, which can be challenging to remove from short signals reliably such as these, potentially impacting inferences about precession and kicks [130].
Finally, the remaining columns in Figs. 11 and 12 report the NRSur7dq4 posteriors and priors for the remnant spin and kick direction parameters. For each parameter (θ χ f , θ v f , ϕ χ f and ϕ v f ), we highlight four interesting events in Fig. 15. These events are picked either because they show the highest JS divergence between posterior and prior, or because the posterior peaks noticeably away from the peak of the prior.
Note that θ χ f = 0 means that the spin is directed along L at t ref = -100 M det . As mentioned above, binaries have a preference for θ χ f = 0 in this frame, and this is reflected in the prior for θ χ f having a strong preference for zero (Fig. 15a). Similarly, Fig. 15b shows that the prior for θ v f has a strong preference for zero or π, which was also motivated above. Finally, to orient the reader, we remind that ϕ χ f = 0 (ϕ v f =0) indicates that the in-plane component of the remnant spin (kick) is coaligned with a vector from the less massive BH to the more massive BH at t ref = -100 M det . In Fig. 15, we note that while the prior for ϕ χ f is broad, the prior for ϕ v f shows a preference for zero. These features highlight why the wave frame at t ref = -100 M det is a suitable frame for discussing remnant spin and kick directions.
In Fig. 12, we find that the JS divergence between the posterior and prior for θ χ f , rises above 0.02 bits for about 91.4% of the events considered, indicating that the data are informative about this parameter for most events. Among the highlighted events in Fig. 15a, GW190630_185205, GW190828_063405 and GW191109_010717 show a stronger preference for θ χ f ∼ 0 than the prior, GW200129_065458 has a slightly bimodal posterior extending to θ χ f ∼ π/6, and GW191109_010717 peaks at θ χ f ∼ π/12, clearly away from the peak of the prior.
Next, for θ v f , we find that JS divergence between the posterior and prior rises above 0.02 bits for about 89.3% of the events considered, indicating that this parameter may be more informative than the kick magnitude v f itself. Among the highlighted events in Fig. 15b, GW170818_022509 prefers a kick directed roughly along L at t ref = -100 M det , while GW200129_065458 prefers the opposite. Here, we note that while GW200129_065458 shows a measurable kick magnitude, GW170818_022509 does not (see Fig. 11). Our findings are broadly consistent with the discussion of the measurability of the kick direction in Refs. [74,140,143], but more work may be needed to understand how to interpret θ v f constraints for events where the kick magnitude is unmeasured.
Finally, the JS divergence between the posterior and prior for ϕ χ f and ϕ v f rises above 0.02 bits for about 6%, and 51% of the events considered, respectively. Therefore, the data are mostly uninformative about ϕ χ f , but can be used to constrain ϕ v f . Among the highlighted events in Fig. 15c and 15d, GW191109_010717 and GW200129_065458 show the strongest constraints; clearly the kick direction is better measured than the kick magnitude.
The third Gravitational-Wave Transient Catalog contains 90 binary coalescence candidates detected by the LIGO-Virgo-KAGRA Collaboration. In this paper, we identify a set of 47 events from the GWTC-3 catalog [5,6,76,77] that falls within the domain of validity (q ≥ 1/6 and M det ≥ 60M ⊙ ) of the NRSur7dq4 waveform model (Sec. III). Within this domain, the NRSur7dq4 model about an order-of-magnitude more accurate than the models used in the the official LVK analysis and includes the full physical effects of precession. We use the Bayesian inference code parallel-bilby to estimate the source properties for these BBH events with the NRSur7dq4 model. We compare source properties inferred using the NRSur7dq4 model and public LVK posteriors obtained using the IMRPhenomXPHM and SEOBNRv4PHM models [5,6,76,77] (Sec. IV). The difference between the resulting posterior samples have been quantified using JS divergence, a common measure of the statistical distance in information content between probability distributions. We find many events for which noticeable differences exist between posteriors obtained with NRSur7dq4, IMRPhenomXPHM, and SEOBNRv4PHM models. Below are the key take-aways from our results:
• While the posteriors are consistent for the majority of events, we find a number of events for which posteriors for NRSur7dq4 are noticeably different from the posteriors for IMRPhenomXPHM/SEOBNRv4PHM. In particular, for ∼ 23% (∼ 55%) of the analyzed events, JS divergence between NRSur7dq4 and IMRPhenomXPHM (SEOBNRv4PHM) exceeds a commonly used threshold of 0.02 bits, indicating nonnegligible differences, for at least one of the following parameters: total mass M , mass ratio q, component masses m 1 , m 2 , spin magnitudes χ 1 , χ 2 , spin tilts θ 1 , θ 2 , the effective inspiral spin χ eff , the spin precession parameter χ p , luminosity distance D L , and inclination angle θ JN (Sec. IV). For many of these cases, multiple parameters exceed this threshold, and for a handful of events, the JS divergence values are above 0.1 (and in a few cases above 0.2) indicating substantial differences. The most interesting GW events are summarized in Sec. IV A.
• Even for the first GW signal, GW150914_095045, we find noticeable differences in χ p and χ 1 measurements between NRSur7dq4, IMRPhenomXPHM and SEOBNRv4PHM. Interestingly, our NRSur7dq4 estimates show better agreement with earlier LVK results obtained with SEOBNRv3 and IMRPhenomPv2.
Fig. 7 summarizes some of these observations.
• For GW191109_010717, NRSur7dq4 shows a stronger preference for negative χ eff at 99.3% credible level, compared to 95.9% for SEOBNRv4PHM and 85.3% for IMRPhenomXPHM (Sec. IV A 6). This is consistent with Ref. [62], where a similar preference was found for the newer SEOBNRv5PHM model. The preference for χ eff < 0 can have important astrophysical implications, as negative χ eff is expected to be more common in dynamically formed binaries than those formed through isolated evolution. However, some caution is warranted as this event suffered from a detector glitch.
• The events showing the most notable differences in the posteriors for the component BH properties are highlighted in Sec. IV A. Furthermore, by comparing the Bayes factors and recovered SNRs between NRSur7dq4 and IMRPhenomXPHM with the same prior settings for all 47 events, we find that there is a mild preference for NRSur7dq4 over IMRPhenomXPHM (Sec. V).
• We find several events where the remnant mass and spin magnitude posteriors are noticeably different between NRSur7dq4 and IMRPhenomXPHM/SEOBNRv4PHM, which can have implications for tests of general relativity (Sec. VI B).
• We provide kick magnitude posteriors for all 47 events, which can be useful for constraining the formation rate of heavy BHs through repeated mergers (Sec. VI C). We find that the kick magnitude is informative for GW191109_010717 and GW200129_065458, with GW200129_065458 showing a preference for a large kick, as noted by Ref. [74]. However, once again, some caution is warranted as both of these events suffered from detector glitches.
• Finally, we also provide posteriors for the remnant spin and kick directions for all 47 events and discuss possible astrophysical applications of these measurements in Sec. VI D.
The differences in the posteriors for NRSur7dq4, IMRPhenomXPHM and SEOBNRv4PHM suggest that waveform systematics are already important for GW data analysis. These differences can become compounded when the posterior samples are used in hierarchical analyses like constraining astrophysical populations or tests of general relativity. Systematic differences can arise from the differences in the modeling approach as well as the physics included; for example, while NRSur7dq4 is trained directly on precessing NR simulations, IMRPhenomXPHM and SEOBNRv4PHM are only informed by nonprecessing simulations.
However, we note that our comparisons are based on posteriors obtained using different Bayesian Inference codes: bilby for IMRPhenomXPHM, RIFT for SEOBNRv4PHM, and parallel-bilby for NRSur7dq4. This may introduce additional systematics in our attempts to make meaningful comparisons. For example, some of the SEOBNRv4PHM posteriors appear undersampled due to an inefficient postprocessing step to include calibration uncertainties in RIFT (Fig. 5). Furthermore, we find that SEOBNRv4PHM posteriors for the skymaps are significantly different for a number of events while NRSur7dq4 and IMRPhenomXPHM match closely (Fig. 9). A more careful study will be necessary to fully disentangle waveform from sampler systematics.
Nonetheless, our results provide further motivation to improve all waveform models. In particular, it is important to extend the region of validity of NR surrogate models to include more unequal mass ratios as well as longer inspirals. As the detector sensitivities improve [152], systematic biases in estimating binary source parameters could limit important applications like BH astrophysics, dark siren cosmology [153], and fundamental tests of general relativity.
Our results are publicly accessible [78]. We additionally provide an application programming interface (API) to access our data programmatically. The API is documented in the GitHub repository for this work [78]. The catalog's website utilized software from the TESS-Atlas project [154].
containing the new live points, shrinking the prior volume that needs to be sampled at subsequent iterations. In contrast, at each iteration of serial nested sampling, only one prior sample is drawn at a time, and a single live point is updated. The differences in how live points are updated affect how bounding ellipses are drawn at each iteration. In practice, we find the parallel variant can prematurely exclude regions of the prior with reasonable posterior support and becomes problematic whenever the number of live-points-per-process becomes too low.
As the number of live-points-per-process increases, the posteriors converge to a unique distribution. Unfortunately, it is not known ahead of time what this number should be. We empirically determine the correct number by randomly selecting five events and systematically varying this value, finding about 16 live-points-per-process is sufficient.
Fig. 16 shows the posterior's dependence as the number of live-points-per-process is varied for one of the representative events GW190727_060333. We consider computations carried out with 640 processes on 1 node (green solid curve) and five nodes (red dashed curve), which give nearly identical posteriors as we would expect based on how parallel-bilby's parallelization is carried out. For an identical setup using 1 node and 128 processes (blue solid curve), we infer a different posterior distribution, In Sec.V, we reanalyzed all 47 events considered in this paper using the IMRPhenomXPHM model, employing the same restricted priors and sampler settings as for the NRSur7dq4 runs (see Sec.II F). These restricted priors, detailed in Sec.II E, are sufficiently broad to encompass the entire posterior distribution. Consequently, comparing posteriors obtained with our setup to the publicly available LVK posteriors [5,6,76,77] using the same IMRPhenomXPHM model serves as a robust test of our parameter estimation (PE) workflow, including sampler settings, PSD computation, data handling, and software versions.
Our primary comparison method involves calculating the Jensen-Shannon divergence (JSD) between various one-dimensional marginalized posteriors. Consistent with the convention used throughout this paper, we interpret JSD values exceeding 0.02 bits as evidence of significant differences between the posteriors from our setup and Jensen-Shannon divergence (JSD) values between the one-dimensional marginalized posteriors for GW150914_095045 for a set of parameters obtained using NRSur7dq4. We compare (shown in green) the public LVK posteriors ("IMRPhenomXPHM") and the posteriors computed using our parameter estimation setup as described in Sections II F, II E, II D, and II C ("IMRPhenomXPHM-RestrictedPrior"). Dashed red lines correspond to a JS divergence of 0.02, indicating significant differences between these posteriors. We find the JDS values (green markers) are well below this threshold, validating the general correctness of our PE setup. For comparison, we show (shown in orange) JSD values for LVK's IMRPhenomXPHM results vs our results obtained with NRSur7dq4. When orange markers are larger than green ones, we can safely ascribe PE differences to waveform systematics as opposed to sampler systematics. Further details on GW150914_095045 are given in Sec. IV A 1.
those published by the LVK. Since the same models and matched parameter estimation (PE) settings are used, we expect most JSD values to remain below this threshold. Fig. 17 presents this comparison for GW150914_095045. The very small JSD values between the LVK posteriors and our own when using the IMRPhenomXPHM model (green) confirm the correctness of our PE setup. In contrast, the JSD values between the NRSur7dq4 run and the LVK posteriors are significantly larger for a subset of parameters, indicating that the differences originate from waveform systematics rather than the PE setup. Similar results hold for the reanalyzed events; for example, over all 47 events, the 95% JSD value for total mass is 0.005 (LVK vs our PE setup when using IMRPhenomXPHM) and 0.04 (NRSur7dq4 posteriors and the LVK posteriors).
An important question arises: which PE results, ours or the LVK's, are "correct"? Unfortunately, this is difficult to resolve. Even within the LVK's analyses, posteri-ors from different models (IMRPhenomXPHM, SEOBNRv4PHM, and combined samples) frequently exhibit JSD values exceeding the 0.02 threshold. At least for when comparing IMRPhenomXPHM and NRSur7dq4, we have controlled for sampler systematics by using the same PE setup as the LVK and verifying the correctness of this setup. Yet any events analyzed in this paper reveal tensions between IMRPhenomXPHM and NRSur7dq4 models when interpreting gravitational wave signals.
Ultimately, model selection using Bayes factors is needed to assess which set of results is more likely to be correct 9 . As reported in Sec. V, the Bayes factors slightly favor the NRSur7dq4 model. However, this preference is not conclusive.
In Sec. V we noted that while there is no outright strong preference for either NRSur7dq4 or IMRPhenomXPHM when considering log Bayes factors or the SNR, the data seems to show a mild preference for NRSur7dq4 over IMRPhenomXPHM. In this appendix, we explore possible explanations for trends and outliers found while performing the model selection.
The NRSur7dq4 model has been trained on a parameter domain defined by q ≥ 1/4 and χ 1 , χ 2 ≤ 0.8, yet throughout this paper, we have used it over the expanded region q ≥ 1/6 and χ 1 , χ 2 ≤ 0.99 (see Sec.II E). We note, however, that IMRPhenomXPHM is also not well calibrated the q < 1/4 or χ 1 > 0.8 region, and wherever comparisons to NR are possible, NRSur7dq4 is more accurate than existing waveform models [67,74,82].
To investigate the possibility that the extrapolated model may produce inaccurate waveforms, we compute the fraction of NRSur7dq4 posteriors that require extrapolation in either the mass ratio or the primary spin magnitude -i.e. q < 1/4 or χ 1 > 0.8, 10 and plot this fraction in Fig. 18 as a function of the differences in the log Bayes factor (left panel) and the median SNR (right panel) between NRSur7dq4 and IMRPhenomXPHM for all 47 events considered. For the events where this fraction is relatively small (< 0.4), we find a mild correlation where NRSur7dq4 is doing better than IMRPhenomXPHM with increasing fractions. For more extreme events (where the fraction is 9 When compared to precessing NR simulations the NRSur7dq4 is more accurate than IMRPhenomXPHM. Hence one could argue that prior model odds (and hence the Bayes factors) should reflect this. We have taken a conservative approach and assume all models are equally likely when setting prior model odds. 10 We do not include extrapolation in χ 2 in this fraction as χ 2 is poorly measured and is therefore prior dominated (see Fig. 2). above 0.4 -these are highlighted in Fig. 18), there is no clear correlation, and with only 4 events its challenging to say anything meaningful. Synthetic NR injection studies could be used to systematically explore model fidelity in these extreme regions of parameter space.
Another possibility is that the NRSur7dq4 model's limited length causes it to miss some of the higher mode content near f low = 20 Hz, resulting in a smaller Bayes Factor or SNR for some events.
As explained in Sec.II F, we use f low = 20 Hz for the overlap integral, and our NRSur7dq4 evaluations return the full length of the surrogate (about 20 orbits). Given this length restriction, the (2, 2) mode for a M det ∼ 60M ⊙ BBH system will start at 20 Hz [67]. This means the higher harmonics will start at multiples of this frequencyfor example, the next most important harmonic, the (3, 3) mode, will start at 30 Hz. And in general, if the (2, 2) mode's starting frequency is f 22 the initial frequency of the (ℓ, m) waveform is m/2 × f 22 Hz.
By using the full length of the surrogate, the (2,2) mode is guaranteed to start below 20 Hz (for M det ≳ 60M ⊙ ), while the higher harmonics with m > 2 may start above 20 Hz. Here, we consider if the missing lower frequency content of certain subdominant modes could impact the Bayes Factor and/or the SNR recovered by the surrogate.
To investigate this possibility, for each event, we randomly choose 5000 samples from the NRSur7dq4 posteriors and compute the starting frequency of the (3, 3) mode for the NRSur7dq4 model in the co-precessing frame. Figure 19 shows the median and 95 percentile values of the (3, 3) mode starting frequency as a function of the differences in the log Bayes factor (left panel) and the median SNR (right panel) between NRSur7dq4 and IMRPhenomXPHM for all 47 events considered. We find no clear correlations suggesting that any missing lowfrequency content in the higher harmonics of NRSur7dq4 may not be the cause for the IMRPhenomXPHM model recovering a higher SNR for GW190517_055101 and a higher Bayes factor for GW190706_222641.
While there could be numerous other causes, another possibility worth noting is that noise fluctuations can also cause the data to prefer one model over another. Determining this will require more observations and/or careful NR synthetic injection studies in realistic noise. As detector sensitivity improves and the number of events increases, comparisons like Figs. 10, 18 and 19 can help determine which waveform models best describe the observed data.
A more comprehensive set of interactive figures is readily found on the NRSurrogate Catalog website[78].
Note that we obtain IMRPhenomXPHM and SEOBNRv4PHM from the public data release associated with the GWTC-3 catalog[5,6,76,77].
The parameter estimation analysis reported in GWTC-1[111] primarily used the LALInference[120] package.
This restricted prior, which is described in Sec. II E, is sufficiently large to contain the full extent of the posterior. Yet the integral appearing in Eq. (6) is carried out over the prior's domain.
We obtain this estimate from the IMRPhenomXPHM runs and add in quadrature to itself to reflect the error estimate for ∆ log e B.
Comparison: in-house vs LVK IMRPhenomXPHM posteriors