GW170104: Observation of a 50-Solar-Mass Binary Black Hole Coalescence at Redshift 0.2

We describe the observation of GW170104, a gravitational-wave signal produced by the coalescence of a pair of stellar-mass black holes. The signal was measured on January 4, 2017 at 10:11:58.6 UTC by the twin advanced detectors of the Laser Interferometer Gravitational-Wave Observatory during their second observing run, with a network signal-to-noise ratio of 13 and a false alarm rate less than 1 in 70,000 years. The inferred component black hole masses are $31.2^{+8.4}_{-6.0}\,M_\odot$ and $19.4^{+5.3}_{-5.9}\,M_\odot$ (at the 90% credible level). The black hole spins are best constrained through measurement of the effective inspiral spin parameter, a mass-weighted combination of the spin components perpendicular to the orbital plane, $\chi_\mathrm{eff} = -0.12^{+0.21}_{-0.30}.$ This result implies that spin configurations with both component spins positively aligned with the orbital angular momentum are disfavored. The source luminosity distance is $880^{+450}_{-390}~\mathrm{Mpc}$ corresponding to a redshift of $z = 0.18^{+0.08}_{-0.07}$. We constrain the magnitude of modifications to the gravitational-wave dispersion relation and perform null tests of general relativity. Assuming that gravitons are dispersed in vacuum like massive particles, we bound the graviton mass to $m_g \le 7.7 \times 10^{-23}~\mathrm{eV}/c^2$. In all cases, we find that GW170104 is consistent with general relativity.


I. INTRODUCTION
The first observing run of the Advanced Laser Interferometer Gravitational-Wave Observatory (LIGO) [1] identified two binary black hole coalescence signals with high statistical significance, GW150914 [2] and GW151226 [3], as well as a less significant candidate LVT151012 [4,5].These discoveries ushered in a new era of observational astronomy, allowing us to investigate the astrophysics of binary black holes and test general relativity (GR) in ways that were previously inaccessible [6,7].We now know that there is a population of binary black holes with component masses 25 M [5,6], and that merger rates are high enough for us to expect more detections [5,8].
Advanced LIGO's second observing run began on November 30, 2016.On January 4, 2017, a gravitationalwave signal was detected with high statistical significance.Figure 1 shows a time-frequency representation of the data from the LIGO Hanford and Livingston detectors, with the signal GW170104 visible as the characteristic chirp of a binary coalescence.Detailed analyses demonstrate that GW170104 arrived at Hanford ∼ 3 ms before Livingston, and originated from the coalescence of two stellar-mass black holes at a luminosity distance of ∼ 3 billion light-years.
GW170104's source is a heavy binary black hole system, with a total mass of ∼ 50 M , suggesting formation in a sub-solar metallicity environment [6].Measurements of the black hole spins show a preference away from being (positively) aligned with the orbital angular momentum, but do not exclude zero spins.This is distinct from the case for GW151226, which had a strong preference for spins with positive projections along the orbital angular momentum [3].The inferred merger rate agrees with previous calculations [5,8], and could potentially be explained by binary black holes forming through isolated binary evolution or dynamical interactions in dense stellar clusters [6].
Gravitational-wave observations of binary black holes are the ideal means to test GR and its alternatives.They provide insight into regimes of strong-field gravity where velocities are relativistic and the spacetime is dynamic.The tests performed with the sources detected in the first observing run showed no evidence of departure from GR's predictions [5,7]; GW170104 provides an opportunity to tighten these constraints.In addition to repeating tests performed in the first observing run, we also test for modifications to the gravitational-wave dispersion relation.Combining measurements from GW170104 with our previous results, we obtain new gravitational-wave constraints on potential deviations from GR.

II. DETECTORS AND DATA QUALITY
The LIGO detectors measure gravitational-wave strain using two dual-recycled Fabry-Perot Michelson interferometers at the Hanford and Livingston observatories [1,10].After the first observing run, both LIGO detectors underwent commissioning to reduce instrumental noise, and to improve duty factor and data quality (see Appendix A in the Supplemental Material [11]).For the Hanford detector, a high-power laser stage was introduced, and as the first step the laser power was increased from 22 W to 30 W to reduce shot noise [10] at high frequencies.For the Livingston detector, the laser power was unchanged, but there was a significant improvement in low-frequency performance mainly due to the mitigation of scattered light noise.
Calibration of the interferometers is performed by inducing test-mass motion using photon pressure from modulated calibration lasers [12,13].The one-sigma calibration uncertainties for strain data in both detectors for the times used in this analysis are better than 5% in amplitude and 3 • in phase over the frequency range 20-1024 Hz.
At the time of GW170104, both LIGO detectors were operating with sensitivity typical of the observing run to date and were in an observation-ready state.Investigations similar to the detection validation procedures for previous events [2,14] found no evidence that instrumental or environmental disturbances contributed to GW170104.

III. SEARCHES
GW170104 was first identified by inspection of lowlatency triggers from Livingston data [15][16][17].An automated notification was not generated as the Hanford detector's calibration state was temporarily set incorrectly in the low-latency system.After it was manually determined that the calibration of both detectors was in a nominal state, an alert with an initial source localization [18,19] was distributed to collaborating astronomers [20] for the purpose of searching for a transient counterpart.Twenty-eight groups of observers covered the parts of the sky localization using ground-and space-based instruments, spanning from γ ray to radio frequencies as well as high-energy neutrinos [21].
Offline analyses are used to determine the significance of candidate events.They benefit from improved calibration and refined data quality information that is unavailable to low-latency analyses [5,14].The second observing run is divided into periods of two-detector cumulative coincident observing time with 5 days of data to measure the false alarm rate of the search at the level where detections can be confidently claimed.Two independently designed matched filter analyses [16,22] used 5.5 days of coincident data collected from January 4, 2017 to January 22, 2017.
These analyses search for binary coalescences over a range of possible masses and by using discrete banks [23][24][25][26][27][28] of waveform templates modeling binaries with component spins aligned or antialigned with the orbital angular momentum [29].The searches can target binary black hole mergers with detector-frame total masses 2 M ≤ M det 100-500 M , and spin magnitudes up to ∼ 0.99.The upper mass boundary of the bank is determined by imposing a lower limit on the duration of the template in the detectors' sensitive frequency band [30].Candidate events must be found in both detectors by the same template within 15 ms [4].This 15-ms window is determined by the 10-ms intersite propagation time plus an allowance for the uncertainty in identified signal arrival times of weak signals.Candidate events are assigned a detection statistic value ranking their relative likelihood of being a gravitational-wave signal: the search uses an improved detection statistic compared to the first observing run [31].The significance of a candidate event is calculated by comparing its detection statistic value to an estimate of the background noise [4,16,17,22].GW170104 was detected with a network matched-filter signal-to-noise ratio (SNR) of 13.At the detection statistic value assigned to GW170104, the false alarm rate is less than 1 in 70,000 years of coincident observing time.
The probability of astrophysical origin P astro for a candidate event is found by comparing the candidate's detection statistic to a model described by the distributions and rates of both background and signal events [8,32,33].The background distribution is analysis dependent, being derived from the background samples used to calculate the false alarm rate.The signal distribution can depend on the mass distribution of the source systems; however, we find that different models of the binary black hole mass distribution (as described in Sec.VI) lead to negligible differences in the resulting value of P astro .At the detection statistic value of GW170104, the background rate in both matched filter analyses is dwarfed by the signal rate, yielding P astro > 1 − (3 × 10 −5 ).
An independent analysis that is not based on matched filtering, but instead looks for generic gravitational-wave bursts [2,34] and selects events where the signal frequency rises over time [35], also identified GW170104.This approach allows for signal deviations from the waveform models used for matched filtering at the cost of a lower significance for signals that are represented by the considered templates.This analysis reports a false alarm rate of ∼ 1 in 20,000 years for GW170104.

IV. SOURCE PROPERTIES
The source parameters are inferred from a coherent Bayesian analysis of the data from both detectors [36,37].As a cross-check, we use two independent model-waveform families.Both are tuned to numericalrelativity simulations of binary black holes with nonprecessing spins, and introduce precession effects through approximate prescriptions.One model includes inspiral spin precession using a single effective spin parameter χ p [38][39][40]; the other includes the generic two-spin inspiral precession dynamics [41][42][43].We refer to these as the effective-precession and full-precession models, respectively [44].The two models yield consistent results.Table I shows selected source parameters for GW170104; unless otherwise noted, we quote the median and symmetric 90% credible interval for inferred quantities.The final mass (or equivalently the energy radiated), final spin and peak luminosity are computed using averages of fits to numerical-relativity results [45][46][47][48][49].The parameter uncertainties include statistical and systematic errors from averaging posterior probability distributions over the two waveform models, as well as calibration uncertainty [37] (and systematic uncertainty in the fit for peak luminosity).Statistical uncertainty dominates the overall uncertainty as a consequence of the moderate SNR.
For binary coalescences, the gravitational-wave frequency evolution is primarily determined by the component masses.For higher mass binaries, merger and ringdown dominate the signal, allowing good measurements of the total mass M = m 1 + m 2 [53][54][55][56][57].For lower mass binaries, like GW151226 [3], the inspiral is more important, providing precision measurements of the chirp mass M = (m 1 m 2 ) 3/5 /M 1/5 [58][59][60][61].The transition between the regimes depends upon the detectors' sensitivity, and GW170104 sits between the two.The inferred component masses are shown in Fig. 2. The form of the twodimensional distribution is guided by the combination of constraints on M and M. The binary was composed of two black holes with masses m 1 = 31.2+8.4 −6.0 M and ).The one-dimensional distributions include the posteriors for the two waveform models, and their average (black).The dashed lines mark the 90% credible interval for the average posterior.The twodimensional plot shows the contours of the 50% and 90% credible regions plotted over a color-coded posterior density function.For comparison, we also show the two-dimensional contours for the previous events [5].
The black hole spins play a subdominant role in the orbital evolution of the binary, and are more difficult to determine.The orientations of the spins evolve due to precession [62,63], and we report results at a point in the inspiral corresponding to a gravitational-wave frequency of 20 Hz [37].The effective inspiral spin parameter χ eff = (m 1 a 1 cos θ LS1 + m 2 a 2 cos θ LS2 )/M is the most important spin combination for setting the properties of the inspiral [64][65][66] and remains important through to merger [67][68][69][70][71]; it is approximately constant throughout the orbital evolution [72,73].Here θ LSi = cos −1 ( L• Ŝi ) is the tilt angle between the spin S i and the orbital angular momentum L, which ranges from 0 • (spin aligned with orbital angular momentum) to 180 • (spin antialigned); i | is the (dimensionless) spin magnitude, which ranges from 0 to 1, and i = 1 for the primary black hole and i = 2 for the secondary.We use the Newtonian angular momentum for L, such that it is normal to the orbital plane; the total orbital angular momentum differs from this because of post-Newtonian corrections.We infer that χ eff = −0.12+0.21 −0.30 .Similarly to GW150914 [5,37,44], χ eff is close to zero with a preference towards being negative: the probability that χ eff < 0 is 0.82.Our measurements therefore disfavor a large total spin positively aligned with the orbital angular momentum, but do not exclude zero spins.
The in-plane components of the spin control the amount of precession of the orbit [62].This may be quantified by the effective precession spin parameter χ p which ranges from 0 (no precession) to 1 (maximal precession) [39].Figure 3 (top) shows the posterior probability density for χ eff and χ p [39].We gain some information on χ eff , excluding large positive values, but, as for previous events [3,5,37], the χ p posterior is dominated by the prior (see Appendix C in the Supplemental Material [11]).No meaningful constraints can be placed on the magnitudes of the in-plane spin components and hence precession.
The inferred component spin magnitudes and orientations are shown in Fig. 3 (bottom).The lack of constraints on the in-plane spin components means that we learn almost nothing about the spin magnitudes.The secondary's spin is less well constrained as the less massive component has a smaller impact on the signal.The probability that the tilt θ LSi is less than 45 • is 0.04 for the primary black hole and 0.08 for the secondary, whereas the prior probability is 0.15 for each.Considering the two spins together, the probability that both tilt angles are less than 90 • is 0.05.Effectively all of the information comes from constraints on χ eff combined with the mass ratio (and our prior of isotropically distributed orientations and uniformly distributed magnitudes) [5].
The source's luminosity distance D L is inferred from the signal amplitude [37,74] 3. Top: Posterior probability density for the effective inspiral and precession spin parameters, χ eff and χp.The one-dimensional distributions show the posteriors for the two waveform models, their average (black), and the prior distributions (green).The dashed lines mark the 90% credible interval for the average posterior.The two-dimensional plot shows the 50% and 90% credible regions plotted over the posterior density function.Bottom: Posterior probabilities for the dimensionless component spins, cS1/(Gm 2 1 ) and cS2/(Gm 2 2 ), relative to the normal of the orbital plane L. The tilt angles are 0 • for spins aligned with the orbital angular momentum and 180 • for spins antialigned.The probabilities are marginalized over the azimuthal angles.The pixels have equal prior probability (1.6 × 10 −3 ); they are spaced linearly in spin magnitudes and the cosine of the tilt angles.Results are given at a gravitational-wave frequency of 20 Hz.
proportional to the distance, but also depends upon the binary's inclination [59,[75][76][77].This degeneracy is a significant source of uncertainty [57,71].The inclination has a bimodal distribution with broad peaks for face-on and face-off orientations (see Fig. 9 if the Supplemental Material [11]).
For GW150914, extensive studies were made to verify the accuracy of the model waveforms for parameter estimation through comparisons with numerical-relativity waveforms [78,79].GW170104 is a similar system to GW150914 and, therefore, it is unlikely that there are any significant biases in our results as a consequence of waveform modeling.The lower SNR of GW170104 makes additional effects not incorporated in the waveform models, such as higher modes [55,80,81], less important.However, if the source is edge on or strongly precessing, there could be significant biases in quantities including M and χ eff [78].Comparison to numericalrelativity simulations of binary black holes with nonprecessing spins [79], including those designed to replicate GW170104, produced results (and residuals) consistent with the model-waveform analysis.

V. WAVEFORM RECONSTRUCTIONS
Consistency of GW170104 with binary black hole waveform models can also be explored through comparisons with a morphology-independent signal model [82].We choose to describe the signal as a superposition of an arbitrary number of Morlet-Gabor wavelets, which models an elliptically polarized, coherent signal in the detector network.Figure 4 plots whitened detector data at the time of GW170104, together with waveforms drawn from the 90% credible region of the posterior distributions of the morphology-independent model and the binary black hole waveform models used to infer the source properties.The signal appears in the two detectors with slightly different amplitudes, and a relative phase shift of approximately 180 • , because of their different spatial orientations [2].The wavelet-and template-based reconstructions differ at early times because the wavelet basis requires high-amplitude, well-localized signal energy to justify the presence of additional wavelets, while the earlier portion of the signal is inherently included in the binary black hole waveform model.
The waveforms reconstructed from the morphologyindependent model are consistent with the characteristic inspiral-merger-ringdown structure.The overlap [58] between the maximum-likelihood waveform of the binary black hole model and the median waveform of the morphology-independent analysis is 87%, consistent with expectations from Monte Carlo analysis of binary black hole signals injected into detector data [34].We also Time-domain detector data (gray), and 90% confidence intervals for waveforms reconstructed from the morphology-independent wavelet analysis (orange) and binary black hole (BBH) models from both waveform families (blue), whitened by each instrument's noise amplitude spectral density.The left ordinate axes are normalized such that the amplitude of the whitened data and the physical strain are equal at 200 Hz.The right ordinate axes are in units of noise standard deviations.The width of the BBH region is dominated by the uncertainty in the astrophysical parameters.
use the morphology-independent analysis to search for residual gravitational-wave energy after subtracting the maximum-likelihood binary black hole signal from the measured strain data.There is an 83% posterior probability in favor of Gaussian noise versus residual coherent gravitational-wave energy which is not described by the waveform model, implying that GW170104's source is a black hole binary.

VI. BINARY BLACK HOLE POPULATIONS AND MERGER RATES
The addition of the first 11 days of coincident observing time in the second observing run, and the detection of GW170104, leads to an improved estimate of the rate density of binary black hole mergers.We adopt two simple representative astrophysical population models: a distribution that is a power law in m 1 and uniform in [83], and a distribution uniform in the logarithm of each of the component masses [5,8].In both cases, we impose m 1 , m 2 ≥ 5 M and M ≤ 100 M [8].Using the results from the first observing run as a prior, we obtain updated rates estimates of R = 103 +110 −63 Gpc −3 yr −1 for the power law, and R = 32 +33 −20 Gpc −3 yr −1 for the

VII ASTROPHYSICAL IMPLICATIONS
uniform-in-log distribution [5].These combine search results from the two offline matched filter analyses, and marginalize over the calibration uncertainty [32].The range for the merger rate that brackets the two distributions, 12-213 Gpc −3 yr −1 , is consistent with the range 9-240 Gpc −3 yr −1 estimated from the first observing run [5,8].Recalculating the rates directly after observing a new event can bias rate estimates, but this bias decreases with increasing event count and is negligible compared to other uncertainties on the intervals.While the median estimates have not changed appreciably, the overall tightening in the credible intervals is consistent with the additional observation time and the increment in the number of events with significant probability of being astrophysical from 3 to 4.
Following the first observing run, we performed a hierarchical analysis using the inferred masses of GW150914, LVT151012 and GW151226 to constrain the binary black hole mass distribution.We assumed the power-law population distribution described above, treating α as a parameter to be estimated, and found α = 2.5 +1.5 −1.6 [5].With the addition of GW170104, α is estimated to be 2.
(see Appendix D in the Supplemental Material [11]); the median is close to the power-law exponent used to infer the (higher) merger rates.
Stars lose mass throughout their lives; to leave a heavy black hole as a remnant they must avoid significant mass loss.Low-metallicity progenitors are believed to have weaker stellar winds and hence diminished mass loss [112].Given the mass of the primary black hole, the progenitors of GW170104 likely formed in a lower metallicity environment Z 0.5Z [6,100,[113][114][115], but low mass loss may also have been possible at higher metallicity if the stars were strongly magnetized [116].
An alternative to the stellar-evolution channels would be binaries of primordial black holes [117][118][119][120]. GW170104's component masses lie in a range for which primordial black holes could contribute significantly to the dark matter content of the Universe, but merger rates in such scenarios are uncertain [118,121].The potential for existing electromagnetic observations to exclude primordial black holes of these masses is an active area of research [119,[122][123][124][125][126][127][128].
Some of the formation models listed above predict merger rates on the order of ∼ 1-10 Gpc −3 yr −1 [85, 87, 92-96, 107, 110, 115].Given that the rate intervals have now tightened and the lower bound (from the uniform-inlog distribution) is ∼ 12 Gpc −3 yr −1 , these channels may be insufficient to explain the full rate, but they could contribute to the total rate if there are multiple channels in operation.Future observations will improve the precision of the rate estimation, its redshift dependence, and our knowledge of the mass distribution, making it easier to constrain binary formation channels.
Gravitational-wave observations provide information about the component spins through measurements of χ eff , and these measurements can potentially be used to distinguish different formation channels.Dynamically assembled binaries (of both stellar and primordial black holes) should have an isotropic distribution of spin tilts, with equal probability for positive and negative χ eff , and a concentration around zero [129].Isolated binary evolution typically predicts moderate ( 45 • ) spin misalignments [130], since the effect of many astrophysical processes, such as mass transfer [131,132] and tides [133,134], is to align spins with the orbital angular momentum.Black hole spins could become misaligned due to supernova explosions or torques during collapse.Large natal kicks are needed to produce negative χ eff by changing the orbital plane [129,130,135].The magnitude of these kicks is currently uncertain [136][137][138][139][140][141] and also influences the merger rate, with high kicks producing lower merger rates in some population-synthesis models [98,100,115,142].For binary neutron stars there is evidence that large tilts may be possible with small kicks [143][144][145][146], and it is not yet understood if similar torques could occur for black holes [138,[147][148][149].The absolute value of χ eff depends on the spin magnitudes.Small values of |χ eff | can arise because the spin magnitudes are low, or because they are misaligned with the orbital angular momentum or each other.The spin magnitudes for binary black holes are currently uncertain, but GW151226 demonstrated that they can be 0.2 [3], and high-mass x-ray binary measurements indicate that the distribution of black hole spins could extend to larger magnitudes [147].For GW170104, we infer χ eff = −0.12+0.21 −0.30 .This includes the possibility of negative χ eff , which would indicate spin-orbit misalignment of at least one component.It also excludes large positive values, and thus could argue against its source forming through chemically homogeneous evolution, since large aligned spins (a i 0.4) would be expected assuming the complete collapse of the progenitor stars [106].The inferred range is consistent with dynamical assembly and isolated binary evolution provided that the positive orbitaligned spin is small (whether due to low spins or mis-alignment) [129,[150][151][152]. Current gravitational-wave measurements cluster around χ eff ∼ 0 (|χ eff | < 0.35 at the 90% credible level for all events; see Fig. 10 in the Supplemental Material [11]) [5].Assuming that binary black hole spins are not typically small ( 0.2), our observations hint towards the astrophysical population favoring a distribution of misaligned spins rather than near orbit-aligned spins [153]; further detections will test if this is the case, and enable us to distinguish different spin magnitude and orientation distributions [154][155][156][157][158][159].

VIII. TESTS OF GENERAL RELATIVITY
To check the consistency of the observed signals with the predictions of GR for binary black holes in quasicircular orbit, we employ a phenomenological approach that probes how gravitational-wave generation or propagation could be modified in an alternative theory of gravity.Testing for these characteristic modifications in the waveform can quantify the degree to which departures from GR can be tolerated given the data.First, we consider the possibility of a modified gravitational-wave dispersion relation, and place bounds on the magnitude of potential deviations from GR.Second, we perform null tests to quantify generic deviations from GR: without assuming a specific alternative theory of gravity, we verify if the detected signal is compatible with GR.For these tests we use the three confident detections (GW150914, GW151226 and GW170104); we do not use the marginal event LVT151012, as its low SNR means that it contributes insignificantly to all the tests [5].

A. Modified dispersion
In GR, gravitational waves are nondispersive.We consider a modified dispersion relation of the form E 2 = p 2 c 2 +Ap α c α , α ≥ 0, that leads to dephasing of the waves relative to the phase evolution in GR.Here E and p are the energy and momentum of gravitational radiation, and A is the amplitude of the dispersion [160,161].Modifications to the dispersion relation can arise in theories that include violations of local Lorentz invariance [162].Lorentz invariance is a cornerstone of modern physics but its violation is expected in certain quantum gravity frameworks [162,163].Several modified theories of gravity predict specific values of α, including massivegraviton theories (α = 0, A > 0) [163], multifractal spacetime [164] (α = 2.5), doubly special relativity [165] (α = 3), and Hořava-Lifshitz [166] and extradimensional [167] theories (α = 4).For our analysis, we assume that the only effect of these alternative theories is to modify the dispersion relation.
To leading order in AE α−2 , the group velocity of gravitational waves is modified as v g /c = 1 + (α − 1)AE α−2 /2 [161]; both superluminal and subluminal propagation velocities are possible, depending on the sign 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 5. 90% credible upper bounds on |A|, the magnitude of dispersion, obtained combining the posteriors of GW170104 with those of GW150914 and GW151226.We use picoelectronvolts as a convenient unit because the corresponding frequency scale is around where GW170104 has greatest amplitude (1 peV h × 250 Hz, where h is the Planck constant).General relativity corresponds to A = 0. Markers filled at the top (bottom) correspond to values of |A| and α for which gravitational waves travel with superluminal (subluminal) speed.
of A and the value of α.A change in the dispersion relation leads to an extra term δΨ(A, α) in the evolution of the gravitational-wave phase [160].We introduce such a term in the effective-precession waveform model [38] to constrain dispersion for various values of α.To this end, we assume flat priors on A. In Fig. 5 we show 90% credible upper bounds on |A| derived from the three confident detections.We do not show results for α = 2 since in this case the modification of the gravitational-wave phase is degenerate with the arrival time of the signal.
There exist constraints on Lorentz invariance violating dispersion relations from other observational sectors (e.g., photon or neutrino observations) for certain values of α, and our results are weaker by several orders of magnitude.However, there are frameworks in which Lorentz invariance is only broken in one sector [168,169], implying that each sector provides complementary information on potential modifications to GR.Our results are the first bounds derived from gravitational-wave observations, and the first tests of superluminal propagation in the gravitational sector.
The result for A > 0 and α = 0 can be reparametrized to derive a lower bound on the graviton Compton wavelength λ g , assuming that gravitons disperse in vacuum in the same way as massive particles [5,7,170].In this case, no violation of Lorentz invariance is assumed.Using a flat prior for the graviton mass, we obtain λ g > 1.5 × 10 13 km, which improves on the bound of 1.0 × 10 13 km from previous gravitational-wave observations [5,7].The combined bound using the three confident detections is λ g > 1.6 × 10 13 km, or for the graviton mass m g ≤ 7.7 × 10 −23 eV/c 2 .

B. Null tests
In the post-Newtonian approximation, the gravitational-wave phase in the Fourier domain is a series expansion in powers of frequency, the expansion coefficients being functions of the source parameters [60,63,171].In the effective-precession model, waveforms from numerical-relativity simulations are also modeled using an expansion of the phase in terms of the Fourier frequency.To verify if the detected signal is consistent with GR, we allow the expansion coefficients to deviate in turn from their nominal GR value and we obtain a posterior distribution for the difference between the measured and GR values [172][173][174][175][176][177].We find no significant deviation from the predictions of GR [5,7].Combined bounds for GW170104 and the two confident detections from the first observing run [5] do not significantly improve the bounds on the waveform phase coefficients.
Finally, we investigate whether the merger-ringdown portion of the detected signal is consistent with the inspiral part [7,178,179].The two parts are divided at 143 Hz, a frequency close to the median inferred (detector-frame) innermost-stable-circular-orbit frequency of the remnant Kerr black hole.For each part, we infer the component masses and spins, and calculate from these the final mass and spin using fits from numerical relativity, as in Sec.IV [45][46][47][48].We then calculate a two-dimensional posterior distribution for the fractional difference between final mass and spin calculated separately from the two parts [7,179].The expected GR value (no difference in the final mass and spin estimates) lies close to the peak of the posterior distribution, well within the 90% credible region.When combined with the posteriors from GW150914, the width of the credible intervals decreases by a factor of ∼ 1.5, providing a better constraint on potential deviations from GR.
In conclusion, in agreement with the predictions of GR, none of the tests we performed indicate a statistically significant departure from the coalescence of Kerr black holes in a quasicircular orbit.

IX. CONCLUSIONS
Advanced LIGO began its second observing run on November 30, 2016, and on January 4, 2017 the LIGO-Hanford and LIGO-Livingston detectors registered a highly significant gravitational-wave signal GW170104 from the coalescence of two stellar-mass black holes.GW170104 joins two other high-significance events [2,3] and a marginal candidate [4] from Advanced LIGO's first observing run [5].This new detection is entirely consistent with the astrophysical rates inferred from the previous run.The source is a heavy binary black hole system, similar to that of GW150914.Spin configurations with both component spins aligned with the orbital angular momentum are disfavored (but not excluded); we do not significantly constrain the component black holes' spin magnitudes.The observing run will continue until mid 2017.Expanding the catalog of binary black holes will provide further insight into their formation and evolution, and allow for tighter constraints on potential modifications to GR.
Further details of the analysis and the results are given in the Supplemental Material [11].Data for this event are available at the LIGO Open Science Center [180].

SUPPLEMENTAL MATERIAL
Appendix A: Noise performance of the detectors Figure 6 shows a comparison of typical strain noise amplitude spectra during the first observing run and early in the second for both of the LIGO detectors [181].For the Hanford detector, shot-noise limited performance was improved above about 500 Hz by increasing the laser power.There are new broad mechanical resonance features (e.g., at ∼ 150 Hz, 320 Hz and 350 Hz) due to increased beam pointing jitter from the laser, as well as the coupling of the jitter to the detector's gravitationalwave channel that is larger than in the Livingston detector.The increase in the noise between 40 Hz and 100 Hz is currently under investigation.For the Livingston detector, significant reduction in the noise between 25 Hz and 100 Hz was achieved mainly by the reduction of the scattered light that re-enters the interferometer.
To date, the network duty factor of the LIGO detectors in the second observing run is about 51% while it was about 43% in the first observing run.The improvement came from better seismic isolation at Hanford, and fine tuning of the control of the optics at Livingston.

Appendix B: Searches
The significance of a candidate event is calculated by comparing its detection statistic value to an estimate of the background noise [4,16,17,22,31].Figure 7 shows the background and candidate events from the offline searches for compact binary coalescences obtained from 5.5 days of coincident data.At the detection statistic value assigned to GW170104, the false alarm rate is less than 1 in 70,000 years of coincident observing time.

Appendix C: Parameter inference
The source properties are estimated by exploring the parameter space with stochastic sampling algorithms [36].Calculating the posterior probability requires the likelihood of the data given a set of parameters, and the parameters' prior probabilities.The likelihood is determined from a noise-weighted inner product between the data and a template waveform [59].Possible calibration error is incorporated using a frequency-dependent spline model for each detector [183].The analysis follows the approach used for previous signals [5,37,44].
A preliminary analysis was performed to provide a medium-latency source localization [184].This analysis used an initial calibration of the data and assumed a (conservative) one-sigma calibration uncertainty of 10% in amplitude and 10 • in phase for both detectors, a reduced-order quadrature model of the effectiveprecession waveform [38,40,185,186] (the most computationally expedient model), and a power spectral den-sity calculated using a parametrized model of the detector noise [82,187].A stretch of 4 s of data, centered on the event, was analysed across a frequency range of 20-1024 Hz.We assumed uninformative prior probabilities [5,37]; technical restrictions of the reduced-order quadrature required us to limit spin magnitudes to < 0.8 and impose cuts on the masses (as measured in the detector frame) such that m det 1,2 ∈ [5.5, 160] M , M det ∈ [12.3, 45.0] M and mass ratio q = m 2 /m 1 ≥ 1/8.The bounds of the mass prior do not affect the posterior, but the spin distributions were truncated.The source position is not strongly coupled to the spin distribution, and so should not have been biased by these limits [18,77].
The final analysis used an updated calibration of the data, with one-sigma uncertainties of 3.8% in amplitude and 2.2 • in phase for Hanford, and 3.8% and 1.9 • for Livingston, and two waveform models, the effectiveprecession model [38,40,186] and the full-precession model [41][42][43].The spin priors were extended up to 0.99.As a consequence of the computational cost of the full-precession model, we approximate the likelihood by marginalising over the time and phase at coalescence as if the waveform contained only the dominant (2, ±2) harmonics [36].This marginalisation is not exact for precessing models, but should not significantly affect signals with binary inclinations that are nearly face on or face off [44].Comparisons with preliminary results from an investigation using the full-precession waveform without marginalisation confirm that this approximation does not impact results.The two waveform models produce broadly consistent parameter estimates, so the overall results are constructed by averaging the two distributions.As a proxy for the theoretical error from waveform modeling, we use the difference between the results from the two approximants [37].A detailed summary of results is given in Table II, and the final sky localization is shown in Fig. 8.
Figure 9 illustrates the distance, and the angle between the total angular momentum and the line of sight θ JN .The latter is approximately constant throughout the inspiral and serves as a proxy for the binary inclination [62,188].The full-precessing model shows a greater preference (after accounting for the prior) for face-on or face-off orientations with θ JN 0 • or 180 • .This leads to the tail of the D L distribution extending to farther distances.There is a preference towards face-on or face-off inclinations over those which are edge on; the probability that | cos θ JN | > 1/ √ Left: Search results from the binary coalescence search described in [16,17,31].The histogram shows the number of candidate events (orange markers) in the 5.5 days of coincident data and the expected background (black lines) as a function of the search detection statistic.The reweighted SNR detection statistic is defined in [31].GW170104 has a larger detection statistic value than all of the background events in this period.At the detection statistic value assigned to GW170104, the search's false alarm rate is less than 1 in 70,000 years of coincident observing time.No other significant candidate events are observed in this time interval.Right: Search results from an independently-implemented analysis [22], where the detection statistic ln L is an approximate log likelihood ratio statistic that is an extension of [182].The two search algorithms give consistent results.
son, the Kullback-Leibler divergence between two equalwidth normal distributions with means one standard deviation apart is 0.5 nat = 0.72 bit.We cannot gain much insight from these spin measurements, but this may become possible by considering the population of binary black holes [193].Figure 10 shows the inferred χ eff distributions for GW170104, GW150914, LVT151012 and GW151226 [5].Only GW151226 has a χ eff (and hence at least one component spin) inconsistent with zero.The others are consistent with positive or negative effective inspiral spin parameters; the probabilities that χ eff > 0 are 0.18, 0.23 and 0.59 for GW170104, GW150914 and LVT151012, respectively.Future analy-sis may reveal if there is evidence for spins being isotropically distributed, preferentially aligned with the orbital angular momentum, or drawn from a mixture of populations [153,155,158,159].
While we learn little about the component spin magnitudes, we can constrain the final black hole spin.The final black hole's dimensionless spin a f is set by the binary's total angular momentum (the orbital angular momentum and the components' spins) minus that radiated away.Figure 11 illustrates the probability distributions for the final mass and spin.To obtain these, we average results from different numerical-relativity calibrated fits for the final mass [46,47] and for the final TABLE II.Parameters describing GW170104.We report the median value with symmetric (equal-tailed) 90% credible interval, and selected 90% credible bounds.Results are given for effective-and full-precession waveform models; the overall results average the posteriors for the two models.The overall results include a proxy for the 90% range of systematic error estimated from the variance between models.More details of the parameters, and the imprint they leave on the signal, are explained in [37].The optimal SNR is the noise-weighted inner product of the waveform template with itself, whereas the matched-filter SNR is the inner product of the template with the data.

Effective precession Full precession Overall
Detector-frame Total mass M det /M 60.0 +5.7 −5.8 59.9 +5.7 −6.9 59.9 +5.7±0.1 −6.5±1.0 Chirp mass M det /M 24.9 +2.5  12.9 +1.7 −1.7 13.0 +1.7±0.0 −1.7±0.0 Matched-filter SNR ρ h|s 13.3 +0.2 spin [45][46][47].The fitting formulae take the component masses and spins as inputs.We evolve the spins forward from the 20 Hz reference frequency to a fiducial near merger frequency, and augment the aligned-spin final spin fits [46,47] with the contribution from the inplane spins before averaging [48].From comparisons with precessing numerical-relativity results, we estimate that systematic errors are negligible compared to the statistical uncertainties.We follow the same approach for fits for the peak luminosity [47,49]; here the systematic errors are larger (up to ∼ 10%) and we include them in the uncertainty estimate [37].We find that a f = 0.64 +0.09 −0.20 .This is comparable to that for previous events [5,37,44], as expected for near equal-mass binaries [194,195], but extends to lower values because of the greater preference for spins with components antialigned with the orbital angular momentum.
The final calibration uncertainty is sufficiently small to not significantly affect results.To check the impact of calibration uncertainty, we repeated the analysis using the effective-precession waveform without marginalising over the calibration.For most parameters the change is negligible.The most significant effect of calibration uncertainty is on sky localization.Excluding calibration uncertainty reduces the 90% credible area by ∼ 2%.Posterior probability densities for the effective inspiral spin χ eff for GW170104, GW150914, LVT151012 and GW151226 [5], together with the prior probability distribution for GW170104.The distribution for GW170104 uses both precessing waveform models, but, for ease of comparison, the others use only the effective-precession model.The prior distributions vary between events, as a consequence of different mass ranges, but the difference is negligible on the scale plotted.population model to the three probable mergers from first observing run [5] and GW170104.We assume that the two-dimensional mass distribution of mergers is the combination of a power law in m 1 and a flat m 2 distribution, with m min = 5 M [196][197][198], and subject to the constraint that M ≤ M max , with M max = 100 M , matching the analysis from the first observing run [5,8].Our sensitivity to these choices for lower and upper cut-off masses is much smaller than the statistical uncertainty in our final estimate of α.The marginal distribution for m 1 is and the parameter α is the power-law slope of the marginal distribution for m 1 at masses m 1 ≤ M max /2.The initial mass function of stars follows a similar powerlaw distribution [83,199], and the mass distribution of companions to massive stars appears to be approximately uniform in the mass ratio q [200][201][202].While the initialfinal mass relation in binary black hole systems is complicated and nonlinear [114,[203][204][205], this simple form provides a sensible starting point for estimating the mass distribution.
Accounting for selection effects and the uncertainty in our estimates of the masses of our four events, and imposing a flat prior on the parameter α [5], we find α = 2.3 +1.3 −1.4 .Our posterior on α appears in Fig. 12.The inferred posterior on the marginal distribution for m 1 appears in Fig. 13; the turnover for m 1 > 50 M is a consequence of our choice of M max = 100 M in Eq. (D2).

Appendix E: Tests of general relativity
The tests of GR use the same algorithm base described in Sec.C [36] for estimation of source parameters, with appropriate modifications to the analytical waveform models [5,7].In the Fourier domain, gravitational waves from a coalescing binary can be described by where ϑ GR are the parameters of the source (e.g., masses and spins) in GR.The tests of GR we perform, except for the inspiral-merger-ringdown consistency test, introduce a dephasing term with an unknown prefactor that captures the magnitude of the deviation from GR.While we modify the phase of the waveform from its GR value, the amplitude is kept unchanged; this is because our analysis is more sensitive to the phase evolution than the amplitude.We use a non-GR template of the form h(f ) = Ã(f ; ϑ GR )e i[Ψ(f ; ϑGR)+δΨ(f ; ϑGR,X modGR )] ,(E2) 12.The posterior distribution for the power-law slope of the massive component of the binary black hole mass distribution, α, described in the main text, using the three probable events from the first observing run [5] and GW170104.We find the median and 90% credible interval are α = 2.3 +1.3 −1.4 .The black line indicates the Salpeter law [83] slope used in the power-law population for estimating binary black hole coalescence rates.where X modGR is a theory-dependent parameter, which is zero in the usual GR templates.To simulate the non-GR waveform, we used the effective-precession model as a base; all the GR and non-GR parameters are assumed unknown and estimated from the data.
With multiple detections it is possible to combine constraints on X modGR to obtain tighter bounds.For a generic parameter ϑ, we compute a combined posterior distribution by combining the individual likelihoods [206].For each event e i we estimate the marginal likelihood density p(e i |ϑ) using a Gaussian kernel density estimator.This gives a simple representation of the likelihood that can be easily manipulated.The combined posterior distribution is computed by multiplying the marginal likelihoods and the chosen prior distribution, p(ϑ|e 1 , . . ., e N ) ∝ p(ϑ) N i=1 p(e i |ϑ). (E3) This is used to compute bounds on ϑ given N detections.We use the three confident detections (GW150914, GW151226 and GW170104) to set combined bounds on potential deviation from GR, except in the case of the inspiral-merger-ringdown consistency test where only GW150914 and GW170104 are used as GW151226 has insufficient SNR from the merger-ringdown to make useful inferences.

Modified dispersion
We have assumed a generic dispersion relation of the form E 2 = p 2 c 2 + Ap α c α , α ≥ 0. To leading order in AE α−2 , the group velocity of gravitational waves is thus modified as v g /c = 1 + (α − 1)AE α−2 /2.The modified dispersion relation results in an extra term to be added to the gravitational-wave phase [160]: Here M det is the redshifted (detector-frame) chirp mass, h is the Planck constant, and D α is a distance measure, where H 0 is the Hubble constant, Ω m and Ω Λ are the matter and dark energy density parameters [52], respectively.
Table III lists the 90% credible upper bounds on the magnitude of A, where the individual and combined bounds for the three confident detections are shown; we see that depending on the value of α and the sign of A, the combined bounds are better than those obtained from GW170104 alone by a factor of ∼ 1-4.5.For all values of α, these bounds are consistent with the uncertainties one might expect for heavy binary black holes using Fisher-matrix estimates on simulated GW150914like signals [161].
For small values of α, it is useful to recast the results in terms of lower bounds on a length scale λ A = hcA 1/(α−2) , which can be thought of as the range (or the screening length) of an effective potential, which is infinite in GR.In Table IV we report the numerical values of these bounds for α < 2. For α = 3, 4, we instead express the bounds as lower limits on the energy scale at which quantum gravity effects might become important, E QG = A −1/(α−2) [207][208][209][210][211].This facilitates the comparison with existing constraints from other sectors, which we show in Table V.
In the subluminal propagation regime, bounds exist from electromagnetic (spectral time lag in gammaray bursts [210]), neutrino (time delay between neutrino and photons from blazar PKS B1424-418 [211]), and gravitational (absence of gravitational Cherenkov radiation [207,209]) sectors.In the superluminal propagation regime, the only existing limits are from the neutrino sector (absence of Bremsstrahlung from electron-positron pairs [208]).The GW170104 constraints are weaker than existing bounds, but are the first constraints on Lorentz violation in the gravitational superluminal-propagation sector.
The posterior distributions for A have long tails, which makes it difficult to accurately calculate 90% limits with a finite number of samples.To quantify this uncertainty on the bounds, for each value of α and sign of A we use Bayesian bootstrapping [212] to generate 1000 instances of the relevant posterior distribution.We find that the 90% credible upper bounds are estimated within an interval whose 90% credible interval width is 20% of the values reported in Table III.
For the (GR) source parameters, to check for the potential impact of errors from waveform modelling, we analysed the data using both the effective-precession model and the full-precession model.However, the fullprecession model was not adapted in time for tests of GR to be completed for this publication.In the first observing run, we performed tests with two different waveform families [5,7]: the effective-precession model [38,40,186], and a nonprecessing waveform model [42,213].We follow the same approach here, and use the same nonprecessing waveform model used for the matched filter search [29].The use of a nonprecessing waveform should give conservative bounds on the potential error from waveform modelling, as some of the differences may come from the failure to include precession effects [44].We find that the numbers so obtained are consistent with the results of the effective-precession model at the tens of percent level.

Parametrized test
The phase evolution of gravitational waves from compact binaries is well understood within GR.The inspiral portion, corresponding to large orbital separation, can be described analytically using the post-Newtonian  Violin plots for the parametrized test, combining posteriors for GW170104 with the two confident detections made in the first observing run, GW150914 and GW151226 [5].This figure has been corrected according to [228], and its 90% bars have been re-positioned to correct a previous plotting error.expansion [63].Modelling the merger dynamics requires the use of numerical-relativity simulations [214][215][216], whereas the post-merger signal is described in black hole perturbation theory as a superposition of damped sinusoids [217][218][219][220]. Accurate analytical waveforms are obtained by tuning the effective-one-body [29,221,222] or phenomenological models [40,223] to numericalrelativity simulations [186,224,225].Given a phase parameter in the phenomenological model whose value in GR is p i , we modify the waveform by introducing new dimensionless parameters δ pi such that p i → p i (1 + δ pi ) [5,7].In the parametrized null test, we freely vary one δ pi at a time (in addition to the other source parameters) to look for deviations from GR.
The bounds on pi obtained from GW170104 are weaker than those from the two confident detections of the first observing run [5].GW151226 had an SNR comparable to GW170104, but it is from a significantly lower mass system [3,5], and hence places better constraints on the inspiral parameters.GW150914 had an SNR twice that of GW170104 (while being of comparable mass), and thus places the best constraints on the late-inspiral and merger-ringdown parameters.Therefore, instead of reporting bounds from GW170104, we provide updated combined bounds, combining the results from the three events.In Fig. 14 we show a violin plot for each of the test parameters.The parameters are plotted (from the left) following the order in which they appear in the post-Newtonian expansion or enter the phenomenological model (the β and α parameters).For all the parameters, the GR solution (δ pi = 0) is contained in the 90% credible interval.

Inspiral-merger-ringdown consistency test
GR is well tested in weak gravitational fields, but fewer tests have been performed in the strong-field regime [163,226,227].It is possible that deviations from the expected behavior of GR only manifest in the most extreme conditions, where spacetime is highly dynamical.The inspiral-merger-ringdown consistency test checks whether the low-frequency, inspiral-dominated portion of the waveform is consistent with the high-frequency, merger-ringdown portion.The two frequency ranges are analysed separately, and the inferred parameters are compared.The test uses the estimated final black hole mass and spin (calculated from the component masses and spins using numerical-relativity fits as detailed in Sec.C) [7,178].If the waveform is compatible with the predictions of GR, we expect that the parameters inferred from the two pieces will be consistent with each other, although the difference will not, in general, be zero because of detector noise.In Fig. 15, we show the posteriors on the fractional difference in the two estimates of the final mass and spin for GW170104 and GW150914, as well as the combined posterior.The difference in the estimates are divided by the mean of the two estimates to produce the fractional parameters that describe potential departures from the GR predictions: ∆a f /ā f for the spin and ∆M f / Mf for the mass [179].These definitions are slightly different from the ones used in our earlier papers [7,178], but serve the same qualitative role [179].Each of the distributions is consistent with the GR value.The posterior for GW170104 is broader, consistent with this event being quieter, and having a lower total mass, which makes it harder to measure the post-inspiral parameters.The width of the 90% credible intervals for the combined posteriors of ∆M f / Mf are smaller than those computed from GW170104 (GW150914) by a factor of ∼ 1.6 (1.3), and the intervals for ∆a f /ā f are improved by a factor of ∼ 1.4 (1.2).
FIG. 4.Time-domain detector data (gray), and 90% confidence intervals for waveforms reconstructed from the morphology-independent wavelet analysis (orange) and binary black hole (BBH) models from both waveform families (blue), whitened by each instrument's noise amplitude spectral density.The left ordinate axes are normalized such that the amplitude of the whitened data and the physical strain are equal at 200 Hz.The right ordinate axes are in units of noise standard deviations.The width of the BBH region is dominated by the uncertainty in the astrophysical parameters.

FIG. 6 .
FIG.6.Comparison of typical noise amplitude spectra of the LIGO detectors in the first observing run (O1) and the early stages of the second observing run (O2).The noise is expressed in terms of equivalent gravitational-wave strain amplitude.Some narrow features are calibration lines(22)(23)(24) Hz for L1, 35-38 Hz for H1, 330 Hz and 1080 Hz for both), suspension fibers' resonances (500 Hz and harmonics) and 60 Hz power line harmonics.

2 FIG. 8 .FIG. 9 .
FIG.8.A Mollweide projection of the posterior probability density for the location of the source in equatorial coordinates (right ascension is measured in hours and declination is measured in degrees).The location broadly follows an annulus corresponding to a time delay of ∼ 3.0 +0.4 −0.5 ms between the Hanford and Livingston observatories.We estimate that the area of the 90% credible region is ∼ 1200 deg 2 .

FIG. 11 .
FIG. 11.Posterior probability density for the final black hole mass M f and spin magnitude a f .The one-dimensional distributions include the posteriors for the two waveform models, and their average (black).The dashed lines mark the 90% credible interval for the average posterior.The twodimensional plot shows the 50% and 90% credible regions plotted over the posterior density function.
FIG.12.The posterior distribution for the power-law slope of the massive component of the binary black hole mass distribution, α, described in the main text, using the three probable events from the first observing run[5] and GW170104.We find the median and 90% credible interval are α = 2.3 +1.3 −1.4 .The black line indicates the Salpeter law[83] slope used in the power-law population for estimating binary black hole coalescence rates.

FIG. 13 .
FIG.13.The posterior probability distribution for the primary component mass m1 of binary black holes inferred from the hierarchical analysis.The black line gives the posterior median as a function of mass, and the dark light grey bands give the 50% and 90% credible intervals.The colored vertical bands give the 50% credible interval from the posterior on m1 from the analyses of (left to right) GW151226, LVT151012, GW170104, and GW150914.The marginal mass distribution is a power law for m1 ≤ 50 M , and turns over for m1 ≥ 50 M due to the constraint on the two-dimensional population distribution that m1 + m2 ≤ 100 M .
FIG.14.Violin plots for the parametrized test, combining posteriors for GW170104 with the two confident detections made in the first observing run, GW150914 and GW151226[5].This figure has been corrected according to[228], and its 90% bars have been re-positioned to correct a previous plotting error.