An improved analysis of GW150914 using a fully spin-precessing waveform model

This paper presents updated estimates of source parameters for GW150914, a binary black-hole coalescence event detected by the Laser Interferometer Gravitational-wave Observatory (LIGO) on September 14, 2015 [1]. Reference presented parameter estimation [2] of the source using a 13-dimensional, phenomenological precessing-spin model (precessing IMRPhenom) and a 11-dimensional nonprecessing effective-one-body (EOB) model calibrated to numerical-relativity simulations, which forces spin alignment (nonprecessing EOBNR). Here we present new results that include a 15-dimensional precessing-spin waveform model (precessing EOBNR) developed within the EOB formalism. We find good agreement with the parameters estimated previously [2], and we quote updated component masses of $35^{+5}_{-3}\mathrm{M}_\odot$ and $30^{+3}_{-4}\mathrm{M}_\odot$ (where errors correspond to 90% symmetric credible intervals). We also present slightly tighter constraints on the dimensionless spin magnitudes of the two black holes, with a primary spin estimate $0.65$ and a secondary spin estimate $0.75$ at 90% probability. Reference [2] estimated the systematic parameter-extraction errors due to waveform-model uncertainty by combining the posterior probability densities of precessing IMRPhenom and nonprecessing EOBNR. Here we find that the two precessing-spin models are in closer agreement, suggesting that these systematic errors are smaller than previously quoted.

] presented parameter estimation of the source using a 13-dimensional, phenomenological precessing-spin model (precessing IMRPhenom) and an 11-dimensional nonprecessing effective-onebody (EOB) model calibrated to numerical-relativity simulations, which forces spin alignment (nonprecessing EOBNR). Here, we present new results that include a 15-dimensional precessingspin waveform model (precessing EOBNR) developed within the EOB formalism. We find good agreement with the parameters estimated previously [Abbott et al. Phys. Rev. Lett. 116, 241102 (2016).], and we quote updated component masses of 35 þ5 −3 M ⊙ and 30 þ3 −4 M ⊙ (where errors correspond to 90% symmetric credible intervals). We also present slightly tighter constraints on the dimensionless spin magnitudes of the two black holes, with a primary spin estimate < 0.65 and a secondary spin estimate < 0.75 at 90% probability. Abbott et al. [Phys. Rev. Lett. 116, 241102 (2016).] estimated the systematic parameter-extraction errors due to waveform-model uncertainty by combining the posterior probability densities of precessing IMRPhenom and nonprecessing EOBNR. Here, we find that the two precessing-spin models are in closer agreement, suggesting that these systematic errors are smaller than previously quoted. DOI: 10.1103/PhysRevX.6.041014 Subject Areas: Astrophysics, Gravitation

I. INTRODUCTION
The detection of the first gravitational-wave (GW) transient, GW150914, by the Laser Interferometer Gravitational-wave Observatory in 2015 [1] marked the beginning of a new kind of astronomy, fundamentally different from electromagnetic or particle astronomy. GW150914 was analyzed using the most accurate signal models available at the time of observation, which were developed under the assumption that general relativity is the correct theory of gravity. The analysis concluded that GW150914 was generated by the coalescence of two black holes (BHs) of rest-frame masses 36 þ5 −4 M ⊙ and 29 þ4 −4 M ⊙ , at a luminosity distance of 410 þ160 −180 Mpc [2]. Throughout this paper, we quote parameter estimates as the median of their posterior probability density, together with the width of the 90% symmetric credible interval.
The GW signal emitted by a binary black hole (BBH) depends on 15 independent parameters: the BH masses and the BH spin vectors (the intrinsic parameters); the inclination and the phase of the observer in the orbital plane, the sky location of the binary (parametrized by two angles, the right ascension and declination), the polarization angle of the GW, the luminosity distance of the binary, and the time of arrival of the GW at the detector (all of which are known as extrinsic parameters). The task of extracting all 15 parameters from interferometric detector data relies on efficient Bayesian inference algorithms and on the availability of accurate theoretical predictions of the GW signal. State-of-the-art numerical-relativity (NR) simulations [3][4][5][6][7][8] can generate very accurate BBH waveforms over a large region of parameter space; however, this region does not yet include (i) binary configurations that have large dimensionless spins (>0.5), extreme mass ratios (<1=3), and many GW cycles (≥40-60), except for a few cases [8][9][10]; nor does it include (ii) systems undergoing significant spin-induced precession of the orbital plane. In practice, parameter estimation requires very many waveform evaluations that span a large region of parameter space, and a purely NR approach is possible if one coarsely discretizes the intrinsic parameters, as has been done for GW150914 [11], or constructs interpolants (surrogates) across NR simulations [12]. However, a continuous sampling of the intrinsic parameter space, even outside regions where NR runs are available, is unfeasible.
The first parameter-estimation study of GW150914 [2] used two such models: an effective-one-body (EOB, Refs. [13,14]) model that restricts spins to be aligned with the orbital angular momentum [15], and a phenomenological model that includes spin-precession effects governed by four effective spin parameters [16]. Here, we present updated parameter estimates using a fully spin-precessing EOB model [17,18], which is parametrized by the full set of BBH properties listed above, including all six BH-spin degrees of freedom, and which reflects additional physical effects described in Sec. II. The inclusion of these effects motivates us to repeat the Bayesian analysis of GW150914 with precessing EOB waveforms. This model was not used in Ref. [2] because it requires costly timedomain integration for each set of BBH parameters; thus, not enough Monte Carlo samples had been collected by the time the study was finalized [19].
The main result of our analysis is that the two precessing models (phenomenological and EOB) are broadly consistent, showing largely overlapping 90% credible intervals for all measured binary parameters, more so than the precessing phenomenological and nonprecessing EOB models compared in Ref. [2]. In that study, the parameter estimates obtained with those two models were combined with equal weights to provide the fiducial values quoted in Ref. [1], and they were differenced to characterize systematic errors due to waveform mismodeling. Because the two precessing models yield closer results, we are now able to report smaller combined credible intervals, as well as smaller estimated systematic errors. Nevertheless, the combined medians cited as fiducial estimates in Ref. [1] change only slightly. In addition, we find that some of the intrinsic parameters that affect BBH evolution, such as the in-plane combination of BH spins that governs precession, are constrained better using the precessing EOB model.
Because precessing-EOB waveforms are so computationally expensive to generate, we cannot match the number of Monte Carlo samples used in Ref. [2]. Thus, we carry out a careful statistical analysis to assess the errors of our summary statistics (posterior medians and credible intervals) due to the finite number of samples. We apply the same analysis to the precessing phenomenological and nonprecessing EOB models, and to their combinations. Although finite-sample errors are a factor of a few larger for the precessing EOB model than for the other two, they remain much smaller than the credible intervals, so none of our conclusions is affected. Last, as a further test of the accuracy and consistency of the two precessing models, we use them to estimate the known parameters of a GW150914-like NR waveform injected into LIGO data. The resulting posteriors are similar to those found for GW150914.
This article is organized as follows. In Sec. II, we discuss the modeling of spin effects in the BBH waveforms used in this paper. In Sec. III, we describe our analysis. We present our results in Sec. IV and our conclusions in Sec. V. Throughout the article, we adopt geometrized units, with G ¼ c ¼ 1.

II. MODELING ORBITAL PRECESSION IN BBH WAVEFORM MODELS
Astrophysical stellar-mass BHs are known to possess significant intrinsic spins, which can engender large effects in the late phase of BBH coalescences: they affect the evolution of orbital frequency, and (if the BH spins are not aligned with the orbital angular momentum) they induce the precession of the orbital plane, modulating the fundamental chirping structure of emitted GWs in a manner dependent on the relative angular geometry of binary and observatory [20]. While measuring BH spins is interesting in its own right, the degree of their alignment and the resulting degree of precession hold precious clues to the astrophysical origin of stellar-mass BBHs [21]: Aligned spins suggest that the two BHs were born from an undisturbed binary star in which both components successively collapsed to BHs; nonaligned spins point to an origin from capture events and three-body interactions in dense stellar environments.
Clearly, the accurate modeling of BH-spin effects is crucial to BBH parameter-estimation studies. Now, even state-of-the-art semianalytical waveform models still rely on a set of approximations that necessarily limit their accuracy. These include finite post-Newtonian (PN) order, calibration to a limited number of NR simulations, rotation to precessing frames, and more. Thus, being able to compare parameter estimates performed with different waveform models, derived under different assumptions and approximations (e.g., in time-vs frequency-domain formulations), becomes desirable to assess the systematic biases due to waveform mismodeling. While observing consistent results does not guarantee the absence of systematic errors (after all, multiple models could be wrong in the same way), the fact that we do not observe inconsistencies does increase our confidence in the models.
Such a comparison was performed in the original parameter-estimation study of GW150914 [2], showing consistency between the precessing phenomenological model and the aligned-spin EOBNR model. This result matched the finding that the BH spins were approximately aligned in GW150914, or that precession effects were too weak to be detected, because of the small number of GW cycles and of the (putative) face-on/face-off presentation of the binary. Nevertheless, it may be argued that the conclusion of consistency remained suspect because only one model in the analysis carried information about the effects of precession; conversely, the estimates of mismodeling systematic errors performed in Ref. [2] were likely increased by the fact that the nonprecessing model would be biased by what little precession may be present in the signal.
The analysis presented in this article, which relies on two precessing-spin waveform families, removes both limitations and sets up a more robust framework to assess systematic biases in future detections where spin effects play a larger role. In the rest of this section, we discuss the features and formulation of the fully precessing EOBNR model. The reader not interested in these technical details (and in the Bayesian-inference setup of Sec. III) may proceed directly to Sec. IV. The precessing EOBNR model (henceforth, "precessing EOBNR") used here can generate inspiral-merger-ringdown (IMR) waveforms for coalescing, quasicircular BH binaries with mass ratio 0.01 ≤ q ≡ m 2 =m 1 ≤ 1, dimensionless BH spin magnitudes 0 ≤ χ 1;2 ≡ jS 1;2 j=m 2 1;2 ≤ 0.99, and arbitrary BH spin orientations [22]. We denote with m 1;2 the masses of the component objects in the binary and with S 1;2 their spin vectors. Note that the model was calibrated only to 38 nonspinning NR simulations that span a smaller portion of the parameter space than defined above, but it was not calibrated to any precessing NR waveform (see below for more details).
The fundamental idea of EOB models consists in mapping the conservative dynamics of a binary to that of a spinning particle that moves in a deformed Kerr spacetime [13,14,[23][24][25][26][27][28], where the magnitude of the deformation is proportional to the mass ratio of the binary. This mapping can be seen as a resummation of PN formulas [29] with the aim of extending their validity to the strong-field regime. As for dissipative effects, EOB models equate the loss of energy to the GW luminosity, which is expressed as a sum of squared amplitudes of the multipolar waveform modes. In the nonprecessing limit, the inspiral-plunge waveform modes are themselves resummations of PN expressions [30][31][32] and are functionals of the orbital dynamics. The ringdown signal is described by a linear superposition of the quasinormal modes [33][34][35] of the remnant BH.
EOB models can be tuned to NR by introducing adjustable parameters at high, unknown PN orders. For the precessing EOB model used in this work, the relevant calibration to NR was carried out in Ref. [15] against 38 NR simulations of nonprecessing-spin systems from Ref. [36], with mass ratios up to 1=8 and spin magnitudes up to almost extremal for equal-mass BBHs and up to 0.5 for unequal-mass BBHs.
Furthermore, information from inspiral, merger, and ringdown waveforms in the test-particle limit were also included in the EOBNR model [37,38]. Prescriptions for the onset and spectrum of ringdown for precessing BBHs were first given in Ref. [17] and significantly improved in Ref. [18].
In the model, the BH spin vectors precess according to when the BH spins are oriented generically, the orbital plane precesses with respect to an inertial observer. The orientation of the orbital plane is tracked by the Newtonian orbital angular momentum L N ≡ μr × _ r, where μ ≡ m 1 m 2 =ðm 1 þ m 2 Þ and r is the relative BH separation. One defines a (noninertial) precessing frame whose z axis is aligned with L N ðtÞ, and whose x and y axes obey the minimum-rotation prescription of Refs. [39,40]. In this frame, the waveform amplitude and phase modulations induced by precession are minimized, as pointed out in several studies [39][40][41][42][43]. Thus, the construction of a precessing EOB waveform consists of the following steps: (i) Compute orbital dynamics numerically, by solving Hamilton's equation for the EOB Hamiltonian, subject to energy loss, up until the lightring (or photon-orbit) crossing; (ii) generate inspiral-plunge waveforms in the precessing frame as if the system were not precessing [15]; (iii) rotate the waveforms to the inertial frame aligned with the direction of the remnant spin; (iv) generate the ringdown signal, and connect it smoothly to the inspiral-plunge signal; (v) rotate the waveforms to the inertial frame of the observer.
A phenomenological precessing-spin IMR model (henceforth, "precessing IMRPhenom") was proposed in Refs. [16,44,45]. These waveforms are generated in the frequency domain by rotating nonprecessing phenomenological waveforms [46] from a precessing frame to the inertial frame of the observer, according to PN formulas that describe precession in terms of Euler angles. The underlying nonprecessing waveforms depend on the BH masses and on the two projections of the spins on the Newtonian angular momentum, with the spin of the BH formed through merger adjusted to also take into account the effect of the in-plane spin components. The influence of the in-plane spin components on the precession is modeled with a single-spin parameter (a function of the two BH spins) and also depends on the initial phase of the binary in the orbital plane. Thus, this model only has four independent parameters to describe the 6 spin degrees of freedom, which is justified by the analysis of dominant spin effects performed in Ref. [44].
While both precessing EOBNR and IMRPhenom models describe spin effects, there are important differences in how they account for precession, which is the main focus of this paper.
(1) In precessing IMRPhenom, the precessing-frame inspiral-plunge waveforms are strictly nonprecessing waveforms, while for precessing EOBNR, some precessional effects are included (such as spin-spin frequency and amplitude modulations) since the orbital dynamics that enters the nonprecessing expressions for the GW modes is fully precessing. (2) The precessing EOBNR merger-ringdown signal is generated in the inertial frame oriented along the total angular momentum of the remnant-the very frame where quasinormal mode frequencies are computed in BH perturbation theory. By contrast, precessing IMRPhenom generates the mergerringdown signal directly in the precessing frame.
(3) The IMRPhenom precessing-frame waveforms contain only the dominant ð2; AE2Þ modes [47], while precessing EOBNR also includes ð2; AE1Þ modes in the precessing frame, although these are not calibrated to NR. (4) In IMRPhenom, the frequency-domain rotation of the GW modes from the precessing frame to the inertial frame is based on approximate formulas (i.e., on the stationary-phase approximation), while precessing EOBNR computes the rotations fully in the time domain, where the formulas are straightforward. (5) In precessing IMRPhenom, the frequency-domain formulas for the Euler angles that parametrize the precession of the orbital plane with respect to a fixed inertial frame involve several approximations: Inplane spin components are orbit averaged; the magnitude of the orbital angular momentum is approximated by its 2PN nonspinning expression; the evolution of frequency is approximated as adiabatic; and the PN formulas that regulate the behavior of the Euler angles at high frequencies are partially resummed. By contrast, precessing EOBNR defines these Euler angles on the basis of the completely general motion of L N ðtÞ; this motion is a direct consequence of the EOB dynamics, and as such, it is sensitive to the full precessional dynamics of the six spin components. A priori, it is not obvious that these approximations will not impact parameter estimation for a generic BBH. However, as far as GW150914 is concerned, Ref. [2] showed broadly consistent results between a precessing and a nonprecessing model; a fortiori, we should expect similar results between two precessing models. Indeed, the GW150914 binary is more probable to be face-off or face-on than edgeon with respect to the line of sight to the detector, and the component masses are almost equal [2]: Both conditions imply that subdominant modes play a minor role.
The nonprecessing models that underlie both approximants were tested against a large catalog of NR simulations [15,46,48], finding a high degree of accuracy in the GW150914 parameter region. However, it is important to bear in mind that these waveform models can differ from NR outside the region in which they were calibrated, and they do not account for all possible physical effects that are relevant to generic BBHs, such us higher-order modes. Finally, neither of the two precessing models has been calibrated to any precessing NR simulation. Thus, we cannot exclude that current precessing models are affected by systematics. References [17,18] compared the precessing EOBNR model to 70 NR runs with mild precession (with mass ratios 1 to 1=5, spin magnitudes up to 0.5, generic spin orientations, and each about 15-20 orbital cycles long) finding sky-location and polarization-averaged overlaps typically above 97% without recalibration.
Since the generation of precessing EOBNR waveforms [at least in the current implementation in the LIGO Algorithm Library (LAL)] is a rather time-consuming process (see [19]), when carrying out parameter-estimation studies with this template family, we introduce a timesaving approximation at the level of the likelihood function. Namely, we marginalize over the arrival time and phase of the signal as if the waveforms contained only ð2; AE2Þ inertial-frame modes since in that case the marginalization can be performed analytically [49]. We have determined that the impact of this approximation is negligible by conducting a partial parameter-estimation study where we do not marginalize over the arrival time and phase. We can understand this physically for GW150914 because in a nearly face-on/face-off binary, the ð2; AE1Þ observer-frame modes are significantly subdominant compared to ð2; AE2Þ modes [50].

III. BAYESIAN INFERENCE ANALYSIS
For each waveform model under consideration, we estimate the posterior probability density [51,52] for the BBH parameters, following the prescriptions of Ref. [2]. To wit, we use the LAL implementation of parallel-tempering Markov chain Monte Carlo and nested sampling [49] to sample the posterior density pðϑjmodel; dataÞ as a function of the parameter vector ϑ: To obtain the likelihood LðdatajϑÞ, we first generate the GW polarizations h þ ðϑ intrinsic Þ and h × ðϑ intrinsic Þ according to the waveform model. We then combine the polarizations into the LIGO detector responses h 1;2 by way of the detector antenna patterns: Finally, we compute the likelihood as the sampling distribution of the residuals [i.e., the detector data d k minus the GW response h k ðϑÞ], under the assumption that these are distributed as Gaussian noise characterized by the power spectral density (PSD) of nearby data [49]: where h·j·i denotes the noise-weighted inner product [53]. The prior probability density pðϑÞ follows the choices of Ref. [2]. In particular, we assume uniform mass priors m 1;2 ∈ ½10; 80M ⊙ , with the constraint m 2 ≤ m 1 , and uniform spin-amplitude priors a 1;2 ¼ jS 1;2 j=m 2 1;2 ∈ ½0; 1, with spin directions distributed uniformly on the two-sphere; and we assume that sources are distributed uniformly in Euclidian volume, with their orbital orientation distributed uniformly on the two-sphere. All the binary parameters that evolve during the inspiral (such as tilt angles between the spins and the orbital angular momentum, θ LS 1;2 ) are defined at a reference GW frequency f ref ¼ 20 Hz. Following Ref. [2], we marginalize over the uncertainty in the calibration of LIGO data [54]. This broadens the posteriors but reduces calibration biases.
To assess whether the data are informative with regard to a source parameter (i.e., where it updates the prior significantly), we perform a Kolmogorov-Smirnov (KS) test. Given an empirical distribution (in our case, the Monte Carlo posterior samples) and a probability distribution (in our case, the prior), the KS test measures the maximum deviation between the two cumulative distributions and associates a p-value to that: For samples generated from the probability distribution against which the test is performed, one expects a p-value around 0.5; pvalues smaller than 0.05 indicate that the samples come from a different probability distribution with a high level of significance-that is, there is only a 5% (or less) chance that the two sets of samples come from the same distribution. The outcomes of our KS tests are only statements about how much the posteriors deviate from the respective priors, but they do not tell us anything about the astrophysical relevance of 90% credible intervals.

IV. RESULTS
The first question that we address is whether parameter estimates derived using the two precessing models (precessing IMRPhenom and precessing EOBNR) are compatible. In particular, we compare posterior medians and 90% credible intervals (the summary statistics used in Ref. [2]) for the parameters tabulated in Table I of Ref. [2], as well as additional spin parameters. The nominal values of the medians and 5% and 95% quantiles for the two models are listed in the "EOBNR" and "IMRPhenom" columns of Table I and Fig. 1. However, it is unclear a priori whether any differences are due to the models themselves or to the imperfect sampling of the posteriors in Markov chain Monte Carlo runs. This is a concern especially for the precessing EOBNR results since the slower speed of EOBNR waveform generation means that shorter chains are available for parameter estimation. To gain trust in our comparisons, we characterize the Monte Carlo error of the medians and quantiles by a bootstrap analysis, as follows.
The Monte Carlo runs for the precessing IMRPhenom model produced an equal-weight posterior sampling bootstrap resamplings, compute summary statistics on each, and measure their variation. However, to improve the representativeness of this analysis given the smaller number of samples in play, we use nine additional equalweight populations, obtained by selecting every (1100 þ i) th sample in the original MCMC run, for i ¼ 1; …; 9. For each of the 1000 Bayesian-bootstrap resamplings, we first choose randomly among the ten equal-weight populations.
Monte Carlo errors are expected to shrink as the inverse square root of the number of samples; this is indeed what we observe, with precessing EOBNR finite-sample errors about ð27 000=2700Þ 1=2 ≈ 3 times larger than for precessing IMRPhenom. Table I and Fig. 1 present the results of this study for several key physical parameters of the source of GW150914. With darker colors, we display the finitesample error estimates on the position of the medians and 5% and 95% quantiles. Lighter colors represent the 90% credible intervals.
Combined estimates.-To account for waveform-mismodeling errors in its fiducial parameter estimates, Ref. [2] cited quantiles for combined posteriors obtained by averaging the posteriors for its two models (in Bayesian terms, this corresponds to assuming that the observed GW signal could have come from either model with equal posterior probability). We repeat the same procedure for the two precessing models, and we show the resulting estimates in the column "Overall" of Table I. Quantiles are more uncertain for the precessing combination because of the larger finite-sampling error of precessing EOBNR. Nevertheless, 90% credible intervals are slightly tighter than cited in Ref. [2]. In the Appendix, we provide a graphical representation of the combined estimates.
Posterior histograms: Masses and spin magnitudes.-We now discuss in some detail the salient features of parameter posteriors. In Figs. 2-6, we show the onedimensional marginalized posteriors for selected pairs of parameters and 90% credible intervals (the dashed lines), as obtained for the two precessing models, as well as the twodimensional probability density plots for the precessing EOBNR model. In Fig. 2, we show the posteriors for the source-frame BH masses m 1;2 : These are measured fairly well, with statistical uncertainties around 10%. In Fig. 3, we show the posteriors for the dimensionless spin magnitudes a 1;2 : The bound on a 1 is about 20% more stringent for precessing EOBNR. This is true even if we account for the larger finite-sampling uncertainty in the precessing EOBNR quantiles (see Table I). The final spin presented in Table I and Fig. 1 was obtained, including the contribution from the in-plane spin components to the final spin [57]; previous publications [1,2] just use the contribution from the aligned components of the spins, which remains sufficient for the final mass computation. Just using the aligned components does not change the precessing EOBNR result but gives a precessing IMRPhenom result of 0.66 þ0.04 −0.06 .
Posterior histograms: Spin directions.- Figure 4 reproduces the disk plot of Ref. [2] for precessing EOBNR. In this plot, the three-dimensional histograms of the dimensionless spin vectors S 1;2 =m 2 1;2 are projected onto a plane perpendicular to the orbital plane; the bins are designed so that each contains the same prior probability mass (i.e., histogramming the prior would result in a uniform shading). It is apparent that the data disfavor large spins aligned or antialigned with the orbital angular momentum, consistently with precessing IMRPhenom results. Because precessing EOBNR favors smaller values of the dimensionless spin magnitudes, the plot is darker towards its We show one-dimensional histograms for precessing EOBNR (red) and precessing IMRPhenom (blue); the dashed vertical lines mark the 90% credible intervals. The two-dimensional density plot shows 50% and 90% credible regions plotted over a colorcoded posterior density function.
FIG. 3. Posterior probability densities for the dimensionless spin magnitudes. (See Fig. 2 for details.) center than its counterpart in Ref. [2]. In agreement with that paper, our analysis does not support strong statements on the orientation of the BH spins with respect to the orbital angular momentum. The spin opening angles (the tilts), defined by cosðθ LS 1;2 Þ ¼ ðS 1;2 ·L N Þ=jS 1;2 j, are distributed broadly. However, the KS test described at the end of Sec. III does indicate some deviation between priors and posteriors, with p-values much smaller than 0.05 for cosðθ LS 1 Þ and cosðθ LS 2 Þ.
Posterior histograms: Effective spin parameters.-In Fig. 5, we show the posteriors of the effective spin combinations χ eff [23,[58][59][60] and χ p [44] defined by where S i⊥ is the component of the spin perpendicular to the orbital angular momentum L N , M is the total observed mass, B 1 ¼ 2 þ 3q=2 and B 2 ¼ 2 þ 3=ð2qÞ, and i ¼ f1; 2g. While χ eff combines the projections of the BH spins onto the orbital angular momentum, χ p depends on their in-plane components and thus relates to precessional effects. Both models have credible intervals for χ eff that contain the value 0 and deviate from the prior significantly. The data provide little information about precession but show a slightly stronger preference for lower values of χ p than expressed by our priors; the deviation is more pronounced for precessing EOBNR. The 90% credible intervals contain the value 0 and extend up to about 0.7 and 0.8 for precessing EOBNR and precessing IMRPhenom, respectively. Thus, precessing EOBNR provides a tighter upper bound.
Posterior histograms: Other spin angles.-To explore other possible differences between the two precessing models, we now consider spin parameters that were not reported in Ref. [2]. In particular, we compute posteriors for θ 12 , the opening angle between the spin vectors, and ϕ 12 , the opening angle between the in-plane projections of the spins. The prior on cos θ 12 is uniform in ½−1; 1, while the prior on ϕ 12 is uniform in ½0; 2π. We show these posteriors in Fig. 6. The θ 12 posteriors deviate appreciably from the prior and are rather similar. By contrast, comparing the opening angle between spin projections onto the orbital plane, ϕ 12 , we find that the precessing EOBNR posterior deviates significantly from the prior (with KS p-value ∈ ½0.0077; 0.075), while the precessing IMRPhenom posterior does not (with KS p-value ∈ ½0.30; 0.60). This result is as it should be since in precessing IMRPhenom binaries with identical projection of the total spin on the orbital plane have identical waveforms. Although the KS p-values suggest that the data provide information about θ 12 and ϕ 12 beyond the prior, we note that the 90% credible intervals for both of these parameters cover approximately 90% of their valid ranges and are indistinguishable for each waveform model.
Spin evolution.-All the source parameters discussed above are measured at a reference frequency of 20 Hz. Exploiting the capability of precessing EOBNR of evolving the BH spin vectors in the time domain, we may address the question of estimating values for the spin parameters at the time of the merger. To do so, we randomly sample 1000 distinct configurations from the precessing EOBNR posteriors, and we evolve them to the maximum EOB orbital frequency (a proxy for the merger in the model). We then produce histograms of the evolved values of χ eff and χ p . We find little if any change between 20 Hz and the merger. Indeed, a KS test suggests that the original and evolved distributions are very consistent, with p-values close to 1.
Comparison with numerical relativity.-The precessing EOBNR waveform model has been tested against NR waveforms using simulations from the SXS catalog [17,18,36]. We can provide a targeted cross-check on the accuracy of precessing EOBNR near GW150914 by performing parameter estimation runs on mock NR signals injected into LIGO data. This test is complementary to an ongoing study of the same nature, which, however, does not employ the precessing EOBNR model used in this paper. We use a new LAL infrastructure [61,62] to inject splineinterpolated and tapered NR waveforms into detector data; spins are defined with respect to the orbital angular momentum at a reference frequency of 20 Hz. All higher harmonics of the GW signal are included up to the l ¼ 8 multipole. At the inclinations used in this study, the impact of modes with l > 2 is small but merits further study; a detailed analysis will be presented in a forthcoming paper. We restrict this investigation to a NR waveform that was computed by the SXS Collaboration using the SpEC [63] code and is available in the public waveform catalogue [64] as SXS:BBH:0308. The intrinsic parameters of the NR waveform q ¼ 0.81, a 1 ¼ 0.34, and a 2 ¼ 0.67 are consistent with the results obtained in Ref. [2], and this waveform agrees well with the detector data.
We can freely choose the angle between the line of sight and the angular momentum of the binary for mock NR signals. Since there is some uncertainty in the binary's inclination, we perform one run near maximum a posteriori probability inclination, ι ¼ 2.856 rad (163.6°), and a second one at the upper bound of the 90% credible interval of the marginal probability density function (PDF) of the inclination, ι ¼ 1.2 rad (68.8°). In Fig. 7, we show the two GW polarizations for the NR waveform and the precessing EOBNR model. The spin magnitudes and the mass ratio were fixed to the NR values. The directions of spins are defined to be the same at the initial time: Tilt angles are 18.8°,149.4°, and the azimuthal angles, defined with respect to the initial separation vector, are 30.9°,38.7°for the primary and secondary BHs, respectively. To quantify the agreement between those waveforms, we compute overlaps averaged over the GW polarization and source sky location, which takes into account the uncertainty in those parameters. The polarization-sky-averaged overlap for MaP inclination is 0.997, and for ι ¼ 1.2 rad (68.8°), overlap is 0.993.
We show results for the run with MaP inclination for the source-frame component masses and effective spins in the left and right panels of Fig. 8. The precessing EOBNR and precessing IMRPhenom models show good agreement in the masses and effective precession spin χ p . The posterior PDFs obtained for the effective aligned spin χ eff are slightly offset. All injected values are found within the 90% credible regions. Results for the inclination chosen at the upper bound of the 90% credible interval of the marginal PDF of the inclination are qualitatively similar to the MaP results, except for the PDF of the effective precession spin, which peaks around χ p ∼ 0.4, noticeably above the injected value but still well inside the 90% credible interval. The inclination is ι ¼ 2.856 rad. The alignment of the precessing EOBNR waveform is obtained from the sky-and polarizationaveraged overlap with the NR waveform. For more details on the alignment procedure, see Ref. [18].

V. CONCLUSIONS
We presented an updated analysis of GW150914 with mass estimates of 35 þ5 −3 M ⊙ and 30 þ3 −4 M ⊙ , and we refined parameter estimates using a generalized, fully precessing waveform model that depends on all 15 independent parameters of a coalescing binary in circular orbit. We find this analysis to be broadly consistent with the results in Ref. [2]. By using the difference between two precessing waveform models as a proxy for systematic errors due to waveform uncertainty, we can compute a more accurate systematic error than what was possible in Ref. [2]. By looking at differences in 5% and 95% quantiles between different waveform models in Fig. 1, one can observe, on average, more consistent values when the two precessing models are compared. In addition, this analysis provides an estimate of the systematic error on precessing spin parameters such as the effective precessing spin χ p and the tilt angles [arccosðŜ 1;2 ·L N Þ], which was not available in Ref. [2]. We have also carefully investigated uncertainties due to the finite numbers of samples used to recreate continuous posterior density functions, and we quantified their effects on quoted credible intervals. As in Ref. [2], the statistical error due to finite signal-to-noise ratio dominates the parameter-estimation error.
While we do recover a tighter limit on the spin magnitude of the most massive member of the binary that created GW150914 (< 0.65 at 90% probability), the recovery of the spin parameters (magnitude and tilt angles) is too broad to hint at whether the black hole binary was formed via stellar binary interactions or dynamical capture [21]. This analysis on the first direct detection by LIGO will be applied to future detections [65], with the aim of getting the most accurate and most precise parameter estimate possible. In particular, binaries that have larger mass asymmetry, that are observed for a longer time, and that are more edge-on than GW150914 will display stronger spin-precession effects. Development, India; the Spanish Ministerio de Economía y Competitividad; the Conselleria d'Economia i Competitivitat and Conselleria d'Educació; Cultura i (left panel) and effective aligned χ eff and effective precessing spins χ p (right panel) for an eventlike NR mock signal close to MaP parameters. In the one-dimensional marginalized distributions, we show the precessing EOBNR (red) and precessing IMRPhenom (blue) probability densities with dashed vertical lines marking 90% credible intervals. The two-dimensional plot shows the contours of the 50% and 90% credible regions of the precessing EOBNR over a color-coded posterior density function. The true parameter values are indicated by a red asterisk or black dot-dashed line. The authors gratefully acknowledge the support of the NSF, STFC, MPS, INFN, CNRS, and the State of Niedersachsen/ Germany for provision of computational resources.

APPENDIX: CREDIBLE INTERVALS FOR THE COMBINED POSTERIORS
To compare directly with the results of Ref. [2], Fig. 9 presents the 90% credible intervals obtained with combined nonprecessing-EOBNR and precessing-IMRPhenom models, and with combined precessing-EOBNR precessing-IMRPhenom models. As in Fig. 1, the darker bands visualize uncertainties due to the finite number of samples, as estimated with the Bayesian bootstrap. FIG. 9. Comparison of parameter estimates obtained by combining the nonprecessing-EOBNR and precessing-IMRPhenom models (as in Ref. [2]; light purple bars at the top) and by combining the precessing-EOBNR and precessing-IMRPhenom models (light green bars at the bottom). We show 90% credible intervals for selected GW150914 source parameters. The darker intervals represent uncertainty estimates for the 5%, 50%, and 95% quantiles (from left to right), as estimated by the Bayesian bootstrap.
Australian National University, Canberra, Australian Capital Territory 0200, Australia