Properties and Astrophysical Implications of the 150 M_⊙ Binary Black Hole Merger GW190521

Observational evidence for IMBHs, usually defined as BHs with mass 10²–10⁵ M_⊙ (see, e.g., Miller & Colbert 2004; van der Marel 2004) has long been sought. IMBHs bridge the gap between stellar BHs and supermassive BHs (SMBHs) and might be the missing link to explain the formation of SMBHs (Volonteri 2010; Greene et al. 2019). IMBHs are predicted to form in the early universe via direct collapse of very massive Population III stars (≳230 M_⊙; e.g., Fryer et al. 2001; Heger et al. 2003; Spera & Mapelli 2017) or through collapse of low angular momentum gas clouds (e.g., Loeb & Rasio 1994; Bromm & Loeb 2003; Begelman et al. 2006; Lodato & Natarajan 2006), perhaps passing through a quasi-star phase (e.g., Begelman 2010; Ball et al. 2011). In stellar clusters, IMBHs are predicted to form via dynamical channels such as runaway collisions (Portegies Zwart & McMillan 2002; Gurkan et al. 2004; Portegies Zwart et al. 2004) and hierarchical mergers of smaller BHs (Miller & Hamilton 2002; O'Leary et al. 2006; Giersz et al. 2015), especially in metal-poor star clusters (Mapelli 2016). However, there is no conclusive observational confirmation of IMBHs in globular clusters and other massive star clusters (Gerssen et al. 2002; Gebhardt et al. 2005; Noyola et al. 2008; Anderson & van der Marel 2010; van der Marel & Anderson 2010; Lützgendorf et al. 2011, 2013; Miller-Jones et al. 2012; Nyland et al. 2012; Strader et al. 2012; Lanzoni et al. 2013; Kızıltan et al. 2017; Perera et al. 2017; Lin et al. 2018; Tremou et al. 2018; Baumgardt et al. 2019; Mann et al. 2019; Zocchi et al. 2019).

Several ultraluminous X-ray sources, defined as those with a total luminosity, assumed isotropic, of ≳10³⁹ erg s⁻¹, have been studied as IMBH candidates (Kaaret et al. 2001; Matsumoto et al. 2001; Strohmayer & Mushotzky 2003; Miller & Colbert 2004; van der Marel 2004; Feng & Soria 2011; Sutton et al. 2012; Mezcua et al. 2015; Mezcua 2017; Wang et al. 2015; Kaaret et al. 2017; Lin et al. 2020), but only a few still support evidence for IMBHs. HLX-1 is possibly the strongest IMBH candidate from electromagnetic data (Farrell et al. 2009; Godet et al. 2009; Servillat et al. 2011; Soria et al. 2012; Webb et al. 2012; Cseh et al. 2015), pointing to an IMBH mass of ∼(0.3–30) × 10⁴ M_⊙.

Several IMBH candidates lie at the centers of dwarf galaxies and are associated with low-luminosity active galactic nuclei (AGNs; Filippenko & Sargent 1989; Filippenko & Ho 2003; Barth et al. 2004, 2005; Greene & Ho 2004, 2007; Seth et al. 2010; Dong et al. 2012; Reines et al. 2013; Baldassare et al. 2015, 2016, 2017; den Brok et al. 2015; Mezcua et al. 2016, 2018). Their estimated masses are close to (or above) the upper edge of the IMBH mass range. The final mass of GW190521 is close to the lower end of the IMBH mass range, in an apparent BH desert covering the mass range of ∼10²–10³ M_⊙. Moreover, this final mass is the first confirmation that IMBHs can form through the merger of two less massive BHs.

The mass of the primary component of GW190521 falls within the range where PI is expected to suppress BH formation. PI develops in a star when the effective production of electron–positron pairs in the stellar core softens the equation of state, removing pressure support (Woosley et al. 2007). This leads to a contraction of the core, raising the internal temperature up to the ignition of oxygen or silicon, and the star becomes unstable. PI is expected to develop in stars with helium-core mass ≳32 M_⊙. For helium cores 32 ≲ M_He/M_⊙ ≲ 64, this instability manifests as pulsational PI (PPI): the star undergoes a number of oscillations that eject material and remove the stellar envelope, bringing the star back to a stable configuration after the resulting mass loss (Barkat et al. 1967; Woosley et al. 2007; Chen et al. 2014; Yoshida et al. 2016). After PPI, the star ends its life with a core-collapse supernova (CCSN) or with direct collapse, leaving a compact object less massive than expected in the absence of PPI (Woosley 2017, 2019). For helium cores 64 ≲ M_He/M_⊙ ≲ 135, PI leads to a complete disruption of the star, leaving no compact object, while for even larger helium cores PI drives a direct collapse to a BH (Fowler & Hoyle 1964; Ober et al. 1983; Bond et al. 1984; Heger et al. 2003; Woosley et al. 2007).

The combined effect of PI and PPI is expected to carve a mass gap in the BH mass function, with lower boundary ∼40–65 M_⊙ and upper boundary ≳120 M_⊙ (Heger et al. 2003; Belczynski et al. 2016; Spera & Mapelli 2017; Woosley 2017, 2019; Giacobbo et al. 2018). The boundaries of the mass gap are highly uncertain because they depend on stellar evolution and on our understanding of CCSNe, PPI supernovae, and PI supernovae (Farmer et al. 2019; Marchant et al. 2019; Stevenson et al. 2019; Mapelli et al. 2020). Several formation channels might populate the mass gap. Below, we will review these channels and attempt to interpret the astrophysical origin of GW190521 in this context and to put constraints on different scenarios.

We describe the detection of GW190521 in a companion paper (Abbott et al. 2020b), where we detail the circumstances of the observation and the detection significance using three different search pipelines. The search results are consistent with a coherent astrophysical signal and inconsistent with an instrumental noise origin for the event. In addition, the strain data are consistent with GW emission from the coalescence of a quasi-circular compact binary system.

In this paper, we begin by assuming that the source is indeed the coalescence of such a binary. In Section 2, we give further details about the Bayesian parameter estimation procedure and the posterior probability distributions that provide estimates of the source's intrinsic and extrinsic parameters. We quantify the evidence for orbital precession due to in-plane component spins and the evidence for the presence of higher-order multipoles beyond the dominant ℓ = 2, m = 2 mode in the data.

In Section 3, we discuss the consistency of the observed signal with the coalescence of a quasi-circular compact binary system. We test the consistency of the residual data, after subtraction of the best-fitting signal, with detector noise, as well as the consistency of the merger and ringdown portions of the signal with expectations from waveform models derived from GR.

In Section 4, we present an estimate of the rate per comoving volume for merger events similar to GW190521. In Section 5, we consider the astrophysical implications of the observation, discussing uncertainties on the PI mass gap and proposed astrophysical channels that might populate the mass gap, including hierarchical merger in stellar cluster environments and stellar mergers. In Section 6, we discuss alternative scenarios for the source of GW190521, including a strongly gravitationally lensed merger, an eccentric BBH, or a primordial BBH; we also exclude a cosmic string cusp or kink or a CCSN as possible sources owing to mismatch with the GW data. Finally, we summarize our observations and consider future prospects in Section 7.

To infer the source properties of GW190521, we analyzed 8 s of data in the LIGO and Virgo detectors around the time of the detection. The data are downsampled from 16,384 to 1024 Hz, as we expect no signal power above several hundred Hz owing to the total mass of GW190521. The parameter estimation analysis is done with two independently developed coherent Bayesian inference pipelines, LALInference (Veitch et al. 2015) and RIFT (Lange et al. 2018; Wysocki et al. 2019), which produce consistent results for the inferred source parameters. Both parameter estimation algorithms assume stationary Gaussian noise characterized by the power spectral density (PSD), which is inferred from the data by the BayesLine algorithm (Littenberg & Cornish 2015). We compute the event likelihood in the frequency domain, integrating over the frequency band 11–512 Hz.

We used three distinct GW signal models of BBH coalescence in our analysis: NRSur7dq4 (NRSur PHM), a surrogate waveform model built by directly interpolating NR solutions (Varma et al. 2019); IMRPhenomPv3HM (Phenom PHM), an inspiral-merger-ringdown waveform model that uses phenomenological frequency domain fits combining post-Newtonian calculations of the GW phase and amplitude (Blanchet et al. 1995, 2005; Damour et al. 2001; Arun et al. 2009; Blanchet 2014) with tuning to NR solutions (Khan et al. 2020); and SEOBNRv4PHM (SEOBNR PHM), an inspiral-merger-ringdown waveform model that is based on the effective-one-body formalism (Buonanno & Damour 1999, 2000) and calibrated to NR (Ossokine et al. 2020).

These three waveform models employ different approaches to reproduce the predictions from analytical relativity and numerical relativity; we expect to see differences in the parameter estimation from these three models as a result of those different approaches, and we can interpret those differences as a form of systematic error associated with the modeling. Note that the effects of the astrophysical environment, such as the presence of gas, on the GW waveform are expected to be negligible (Fedrow et al. 2017) in the late stages of inspiral, merger, and ringdown that we observe.

NRSur PHM is constructed based on NR simulations with component spins that are not constrained to be aligned with the orbital axis, thus including the effects of spin–orbit precession. The model covers dimensionless spin magnitudes χ_i ≤ 0.8 and mass ratios q = m₂/m₁ ≥ 1/4. It includes all $(l,| m| )$-multipoles of the gravitational radiation (Blanchet et al. 1996, 2008; Kidder 2008; Mishra et al. 2016) up to and including l = 4. In the training parameter space, the waveform model has shown excellent agreement with NR simulations, with mismatches comparable to the numerical errors associated with the NR simulation. The model continues to agree with NR when extrapolating to mass ratios q ≥ 1/6 (Varma et al. 2019). NRSur PHM is directly trained with NR simulations and therefore only models the last ∼20 orbits of the inspiral, which is adequate for GW190521 because the signal is in the measurement band of the detectors for fewer cycles.

Phenom PHM (Khan et al. 2019) is an approximate higher-multipole aligned-spin waveform model that maps the subdominant radiative moments $(l,| m| )=(2,1),(3,2),(3,3),(4,3),(4,4)$ to the dominant (2, 2) mode (London et al. 2018). Multipoles are defined in the co-precessing frame, where the binary approximates a system with aligned spins (Schmidt et al. 2011; Pekowsky et al. 2013), and are then transformed by a time-dependent rotation to model the harmonic modes of a precessing binary in the inertial frame (Schmidt et al. 2012; Hannam et al. 2014), using a double-spin model of spin–orbit precession during the inspiral (Chatziioannou et al. 2017). After this precession "twisting," all l ≤ 4 modes will be nonzero in the inertial frame. The nonprecessing, dominant multipole of the radiation is tuned to spin-aligned NR simulations in the parameter space of spin magnitudes χ_i ≤ 0.85 and mass ratios q ≥ 1/18 (Husa et al. 2016; Khan et al. 2016). The subdominant multipoles and precessional effects in Phenom PHM, however, have not been calibrated to NR, and the model does not include spherical−spheroidal mode-mixing effects, which can significantly impact some of the higher multipoles (Kelly & Baker 2013).

SEOBNR PHM is based on the dynamics of spinning, nonprecessing BBHs in the effective-one-body formalism, calibrated to NR simulations and results from BH perturbation theory (Bohé et al. 2017). The model includes the nonprecessing $(l,| m| )=(2,1),(3,3),(4,4),(5,5)$ multipoles in addition to the dominant (2, 2) multipole mode. The individual modes are calibrated to waveforms from NR and BH perturbation theory (Cotesta et al. 2018), covering the parameter space of mass ratios q ≥ 1/10 and effective inspiral spin parameter χ_eff ∈ [−0.7, 0.85]. Effects from the precessing orbital plane are modeled through a suitable rotation of the nonprecessing inspiral-plunge multipoles from the co-precessing frame to the inertial frame (without recalibration to NR), with a direct attachment of merger-ringdown modes in the co-precessing frame (Babak et al. 2017; Ossokine et al. 2020). After the precession twisting, all l ≤ 5 modes will be nonzero in the inertial frame.

Although not directly fitted to NR simulations of precessing binaries, both Phenom PHM and SEOBNR PHM have been validated through a comparison with a large set of such NR waveforms (Khan et al. 2020; Ossokine et al. 2020). The three models described above are tuned to different NR solutions. Comparisons between different NR codes find agreement on the level of accuracy of the individual codes (Hannam et al. 2009; Hinder et al. 2014; Lovelace et al. 2016). The agreement between different NR codes is sufficiently good (Pürrer & Haster 2020) to avoid systematic biases at the S/N of GW190521. However, The three models are constructed in sufficiently different ways that it is useful to compare the results of parameter estimation from each of them.

The NRSur PHM waveform model is most faithful to NR simulations in the parameter range relevant for GW190521 (Varma et al. 2019). Therefore, the inferred source parameters quoted in this paper and in Abbott et al. (2020b) were obtained with the NRSur PHM model unless otherwise noted. In this section we also present results from the Phenom PHM and SEOBNR PHM models to check for systematic differences between waveform models. As shown below, differences in results between waveform models are of the same order as or smaller than the statistical error and do not impact the astrophysical interpretation of GW190521 discussed in this paper.

We choose priors that are uniform on the component masses in the detector frame from [30, 200] M_⊙. We further restrict the mass priors such that the total mass must be greater than 200 M_⊙ and the chirp mass to be between 70 and 150 M_⊙, both in the detector frame. In all cases, we verify that the posterior distributions do not have support at the boundaries of the priors. The distance prior scales as ${D}_{{\rm{L}}}^{2}$ (i.e., differential in Euclidean volume) up to a maximum of 10 Gpc. We have checked that a prior that is differential in comoving volume (with Planck 2015 cosmology) makes negligible (<1%) difference in the posterior medians for all source parameters. For the BH spins, we adopt a uniform prior for the magnitude in the dimensionless parameters χ_i ∈ [0, 0.99] and their orientation angles chosen to be uniform on the surface of the unit sphere. We adopt a uniform prior in the cosine of the inclination angle between the binary angular momentum and the line of sight, θ_JN. The prior on the sky location (R.A. and decl.) is chosen to be uniform on the surface of the unit sphere.

The parameters of the remnant BH formed after the merger, its mass M_f, dimensionless spin χ_f, and recoil kick velocity v_f, are inferred by applying fits calibrated to NR to the posterior distributions of the binary's initial masses and spins. For the posteriors from the Phenom PHM and SEOBNR PHM waveform models, we used the same M_f and χ_f fits that are implemented internally in these models: for Phenom PHM, the fits from Husa et al. (2016) with corrections for precession from Bohé et al. (2016), and for SEOBNR PHM, the fits from Ossokine et al. (2020) and from Hofmann et al. (2016) applied to the spins evolved using the waveform model's dynamics as described in Ossokine et al. (2020). For the NRSur PHM posterior, we applied the related surrogate remnant fit of Varma et al. (2019) for M_f, χ_f, and v_f. Applying the M_f and χ_f surrogate fits to the posteriors from the Phenom PHM and SEOBNR PHM waveform models, and using the average of the fits from Healy & Lousto (2017), Hofmann et al. (2016), and Jiménez-Forteza et al. (2017) after applying corrections for precession (Johnson-McDaniel et al. 2016; Abbott et al. 2017b), both yield consistent results. The surrogate v_f was only tested for NRSur PHM (Varma et al. 2020); therefore, we do not apply it to Phenom PHM and SEOBNR PHM. The peak luminosity is also inferred using fits calibrated to NR (Healy & Lousto 2017; Keitel et al. 2017), while the energy radiated in the merger is given by M − M_f. The key analysis elements described above, including parameter estimation sampling algorithms, PSD estimates, and waveform models, all potentially introduce systematic uncertainties. Different choices for these elements can affect the results but in most cases these changes are significantly smaller than the statistical uncertainties. Below, we highlight the more significant differences in the results associated with waveform models.

In Table 1 we summarize the source properties of GW190521. Results are quoted as the median and symmetric 90% credible interval of the marginalized posterior distributions for each parameter, and for each of the three GW signal models. The measurements are marginalized over uncertainty in the data calibration. In the rest of this paper we quote source properties derived using NRSur PHM, unless explicitly stated otherwise.

Table 1.  Source Properties for GW190521: Median Values with 90% Credible Intervals That Include Statistical Errors

Waveform Model	NRSur PHM	Phenom PHM	SEOBNR PHM
Primary BH mass m1 (M⊙)	${85}_{-14}^{+21}$	${90}_{-16}^{+23}$	${99}_{-19}^{+42}$
Secondary BH mass m2 (M⊙)	${66}_{-18}^{+17}$	${65}_{-18}^{+16}$	${71}_{-28}^{+21}$
Total BBH mass M (M⊙)	${150}_{-17}^{+29}$	${154}_{-16}^{+25}$	${170}_{-23}^{+36}$
Binary chirp mass ${ \mathcal M }$ (M⊙)	${64}_{-8}^{+13}$	${65}_{-7}^{+11}$	${71}_{-10}^{+15}$
Mass ratio q = m2/m1	${0.79}_{-0.29}^{+0.19}$	${0.73}_{-0.29}^{+0.24}$	${0.74}_{-0.42}^{+0.23}$
Primary BH spin χ1	${0.69}_{-0.62}^{+0.27}$	${0.65}_{-0.57}^{+0.32}$	${0.80}_{-0.58}^{+0.18}$
Secondary BH spin χ2	${0.73}_{-0.64}^{+0.24}$	${0.53}_{-0.48}^{+0.42}$	${0.54}_{-0.48}^{+0.41}$
Primary BH spin tilt angle ${\theta }_{{\mathrm{LS}}_{1}}\,(\deg )$	${81}_{-53}^{+64}$	${80}_{-49}^{+64}$	${81}_{-45}^{+49}$
Secondary BH spin tilt angle ${\theta }_{{\mathrm{LS}}_{2}}\,(\deg )$	${85}_{-55}^{+57}$	${88}_{-58}^{+63}$	${93}_{-60}^{+61}$
Effective inspiral spin parameter χeff	${0.08}_{-0.36}^{+0.27}$	${0.06}_{-0.39}^{+0.31}$	${0.06}_{-0.35}^{+0.34}$
Effective precession spin parameter χp	${0.68}_{-0.37}^{+0.25}$	${0.60}_{-0.44}^{+0.33}$	${0.74}_{-0.40}^{+0.21}$
Remnant BH mass Mf (M⊙)	${142}_{-16}^{+28}$	${147}_{-15}^{+23}$	${162}_{-22}^{+35}$
Remnant BH spin χf	${0.72}_{-0.12}^{+0.09}$	${0.72}_{-0.15}^{+0.11}$	${0.74}_{-0.14}^{+0.12}$
Radiated energy Erad (M⊙ c2)	${7.6}_{-1.9}^{+2.2}$	${7.2}_{-2.2}^{+2.7}$	${7.8}_{-2.3}^{+2.8}$
Peak Luminosity ℓpeak (erg s−1)	${3.7}_{-0.9}^{+0.7}$ × 1056	${3.5}_{-1.1}^{+0.7}\times {10}^{56}$	${3.5}_{-1.4}^{+0.8}\times {10}^{56}$
Luminosity distance DL (Gpc)	${5.3}_{-2.6}^{+2.4}$	${4.6}_{-1.6}^{+1.6}$	${4.0}_{-1.8}^{+2.0}$
Source redshift z	${0.82}_{-0.34}^{+0.28}$	${0.73}_{-0.22}^{+0.20}$	${0.64}_{-0.26}^{+0.25}$
Sky localization ${\rm{\Delta }}{\rm{\Omega }}\,({\deg }^{2})$	774	862	1069

Download table as:  ASCIITypeset image

Masses.—The estimated mass posterior distributions are shown in the left panel of Figure 1 for the three GW signal models. The primary BH mass of GW190521 is m₁ = ${85}_{-14}^{+21}$ M_⊙, making it the highest-mass component BH known to date in GW astronomy. The mass of the secondary BH is inferred to be m₂ = ${66}_{-18}^{+17}$ M_⊙. The primary BH of GW190521 is more massive (median value) than any remnant BH reported in GWTC-1 except for GW170729 (Abbott et al. 2019i); the secondary BH of GW190521 is also more massive than any primary BH in GWTC-1.

Zoom In Zoom Out Reset image size

These source-frame masses have been redshift corrected, as discussed above, using a value of the Hubble parameter H₀ = 67.9 from Planck 2015. However, recent measurements of H₀ using nearby Cepheid distance standards obtain a precise value of H₀ = 74.03 ± 1.42 km s⁻¹ Mpc⁻¹ (Riess et al. 2019), 9% higher than the Planck value. Using the latter value along with the other cosmological parameters from Planck 2015 increases the median value of the redshift by 7% and reduces the estimated source-frame masses by 3%. These shifts are significantly smaller than statistical or other systematic uncertainties, including those affecting the astrophysical interpretation discussed throughout this paper.

While the low-mass cutoff of the PI mass gap is uncertain (see Section 5.1), the primary BH of GW190521 offers strong evidence for the existence of BHs in the mass gap. If the PI gap begins at 50 M_⊙ (65 M_⊙), we find that the primary BH has only a < 0.1% (0.3%) probability of being below the mass gap, while the secondary BH has 6.6% (46.2%) probability of also being below the mass gap.

The SEOBNR PHM model supports a higher primary mass and more asymmetric mass ratio for GW190521: within 90% credible intervals, m₁ and m₂ can be as high as 141 M_⊙ and 92 M_⊙ respectively, while support for the mass ratio extends down to q ∼ 0.32. While the upper limit of the PI mass gap remains uncertain, adopting 120 M_⊙ as the high-mass end of the gap, we find the probability that the primary BH of GW190521 is beyond the gap of 12% when using the SEOBNR PHM model. The corresponding probabilities using the NRSur PHM and Phenom PHM models are 0.9% and 2.3%, respectively.

The probability that at least one of the BHs in GW190521 is in the range 65–120 M_⊙ is 99.0%, using the NRSur PHM model. The corresponding probabilities using the SEOBNR PHM and Phenom PHM models are 90.2% and 98.0%, respectively.

We measure the total binary mass of GW190521 to be M = ${150}_{-17}^{+29}$ M_⊙ making it the highest-mass binary observed via GWs to date. The binary chirp mass is ${ \mathcal M }={64}_{-8}^{+13}$ M_⊙, a factor of ∼2 higher than the first BBH detection, GW150914 (Abbott et al. 2016a, 2019i). GW190521 is consistent with a nearly equal mass binary with mass ratio q = m₂/m₁ = ${0.79}_{-0.29}^{+0.19}$ (90% credible interval).

In the detector frame, the measured masses are ${m}_{1}^{det}\,={152}_{-19}^{+32}$ M_⊙, ${m}_{2}^{det}={120}_{-32}^{+21}$ M_⊙, M^det = ${273}_{-27}^{+26}$ M_⊙, and ${{ \mathcal M }}^{\det }\,={117}_{-14}^{+12}$ M_⊙, using the NRSur PHM model. These results are very nearly the same for all three models.

Spins.—Due to its high total mass, GW190521 is the shortest-duration signal (approximately 0.1 s) recorded so far in the LIGO and Virgo detectors. With only around four cycles (two orbits) in the frequency band 30–80 Hz (Abbott et al. 2020b), information about spin evolution during the coalescence is limited. Still, analyses of GW190521 indicate that GW signal models including effects of spin–orbit precession are mildly preferred over those that omit such effects (i.e., allow only spins aligned with the orbital axis), with a ${\mathrm{log}}_{10}$-Bayes factor of 1.06 ± 0.06 for the NRSur PHM model allowing generic BH spins versus limiting the effects of spin to the aligned components.

In the disk plots of Figure 2, we show constraints on the spins of the component BHs of GW190521 in terms of their dimensionless magnitudes χ₁ and χ₂ and polar angles (tilts) with respect to the orbital angular momentum, ${\theta }_{{\mathrm{LS}}_{1}}$ and ${\theta }_{{\mathrm{LS}}_{2}}$, defined at a fiducial GW frequency of 11 Hz. Median values from all three waveform models suggest in-plane spin components with high spin magnitudes for both the BHs. Within the 90% credible intervals given in Table 1, however, the constraints on the dimensionless BH spin magnitudes remain uninformative. For our preferred model NRSur PHM, the 90% bounds on spin magnitude extend from χ_1,2 ∼ 0.1 to 0.9. The constraints on the tilt angles of these spins are also relatively broad.

Zoom In Zoom Out Reset image size

As for past GW observations, we present inferences on the spins of GW190521 using the parameters χ_eff and χ_p constructed from the mass and spin of the binary components. Here, ${\chi }_{\mathrm{eff}}$ = $({m}_{1}{\chi }_{1}\cos {\theta }_{1}+{m}_{2}{\chi }_{2}\cos {\theta }_{2})/({m}_{1}+{m}_{2})\in [-1,1]$ is the effective inspiral spin parameter (Damour 2001; Ajith et al. 2011), which measures the mass-weighted net spin aligned with the orbital angular momentum axis $\hat{{\boldsymbol{L}}}$ and remains approximately constant throughout the inspiral (Racine 2008). The effective precession spin parameter χ_p ∈ [0, 1] (Schmidt et al. 2015) measures the spin components in the plane of the orbit and therefore the strength of the spin–orbit precession in the binary. The inferred posterior distributions of χ_eff and χ_p are shown in the right panel of Figure 1 for the three GW signal models. Our priors are uniform in the component spins of the binaries, which results in nontrivial priors for the effective spin parameters, shown as the black distribution in the figure. For all three waveform models the posterior distribution of χ_eff is peaked close to 0, similarly to the prior, while the χ_p distribution is shifted toward higher values. While the bulk of the posterior on χ_p suggests an in-plane component of the spins, which contributes to spin-induced precession, the broad distribution prevents a more conclusive finding.

Mass and Spin.—The merger of GW190521 resulted in a final (remnant) BH of mass M_f = ${142}_{-16}^{+28}$ M_⊙ (see Figure 3). The inferred mass of the remnant BH provides observational evidence for an IMBH of ≳100 M_⊙. The remnant BH mass is ${7.6}_{-1.9}^{+2.2}$ M_⊙ less than the sum of the component BH masses; the equivalent energy was released as GWs during coalescence, making GW190521 the most energetic GW event recorded to date. We find a peak luminosity close to merger ℓ_peak = ${3.7}_{-0.9}^{+0.7}$ × 10⁵⁶ erg s⁻¹.

Zoom In Zoom Out Reset image size

The remnant BH of GW190521 has a dimensionless spin parameter χ_f = ${0.72}_{-0.12}^{+0.09}$. Within the 90% credible interval this is consistent with the BBH merger remnant spins reported in GWTC-1 (Abbott et al. 2019i). The predictions for remnant BH parameters from the inferred values of the component masses and spins agree with analyses of GW190521 that target the ringdown portion of the signal, directly measuring M_f and χ_f without assuming a quasi-circular BBH, described in Section 3.2.

Recoil Velocity.—In generic BBH mergers, the radiation of linear momentum through beamed GW emission imparts a recoil velocity, or kick, to the remnant BH (Fitchett 1983; Favata et al. 2004), of magnitude up to ∼5000 km s⁻¹ for precessing systems (Favata et al. 2004; Lousto & Zlochower 2011). As large in-plane spin components are not ruled out for GW190521, it is a potential candidate for a large kick. Figure 4 shows the prior and posterior distributions for the kick magnitude, with respect to the center of mass of the progenitor binary, computed using NRSur PHM and applying the related remnant surrogate model (Varma et al. 2019) to the component masses and spins. Although the kick velocity remains unconstrained, we find that the posterior is weighted toward higher values relative to the prior, with support for values exceeding fiducial escape velocities for globular clusters, the half-mass radius of nuclear star clusters (Antonini & Rasio 2016), giant elliptical galaxies (Fitchett 1983; Merritt et al. 2004; Campanelli et al. 2007a), and Milky Way–like galaxies (Monari et al. 2018).

Zoom In Zoom Out Reset image size

Sky Position.—The initial sky map of GW190521, computed using the low-latency pipeline BAYESTAR (Singer & Price 2016), was publicly released with the initial public circular on 2019 May 21 03:08:32 UTC (Abbott et al. 2019b). The source was localized within a 90% credible area of 1163 deg² (see blue curve in Figure 5). An updated sky map was released the same day at 13:32:27 UTC (Abbott et al. 2019f), using a low-latency analysis from LALInference employing the SEOBNRv4_ROM model (Bohé et al. 2017), giving a 90% credible localization within 765 deg² (orange curve in Figure 5). Here we report our latest constraints on the sky position of GW190521 found using LALInference with the NRSur PHM model (green curve in Figure 5). The source is now localized within 774 deg².

Zoom In Zoom Out Reset image size

Distance and Inclination.—The inferred luminosity distance D_L reported at low latency (Abbott et al. 2019b, 2019f) made use of waveform models that did not include higher-order multipoles. As discussed in Section 2.5, the luminosity distance and redshift inferred using the waveform models described in this paper that incorporate higher-order multipoles result in larger values than those low-latency estimates. In Figure 6, we show the posteriors on the luminosity distance D_L and the inclination angle θ_JN of GW190521 for the three waveform models. Here, the inclination angle is between the total angular momentum of the source and the observer's line of sight. One can note in the figure the typical correlation between D_L and θ_JN and near degeneracy between systems with the angular momentum vector pointed toward or away from the observer. GW emission is strongest along the orbital angular momentum direction, so face-on sources at larger distances produce similar signals to edge-on sources closer by. Both the NRSur PHM and Phenom PHM models suggest that the total angular momentum of the source is roughly aligned with the observer's line of sight θ_JN ∼ 0° or 180°, i.e., the orbital plane is close to face-on to the line of sight, while the SEOBNR PHM model also supports an orbital plane that is closer to edge-on (θ_JN ∼ 90°). In accordance with the covariance between the distance and inclination, the NRSur PHM model places GW190521 at D_L = ${5.3}_{-2.6}^{+2.4}$ Gpc (z ≃ ${0.82}_{-0.34}^{+0.28}$), while the SEOBNR PHM model suggests D_L = ${4.0}_{-1.8}^{+2.0}$ Gpc. The masses of the component BH and remnant BH scale inversely with (1 + z). In accordance, as reported in Table 1, the SEOBNR PHM model supports higher masses than the other two models. A binary close to face-on, as favored by our preferred NRSur PHM model, makes it difficult to directly observe the effects of spin–orbit precession; thus, the evidence for precession in GW190521, reported in the Spins subsection of Section 2.2 above, is weak.

Zoom In Zoom Out Reset image size

The GWs emitted during a BBH coalescence can be decomposed as spherical harmonic multipoles ℓ and $| m| \leqslant {\ell }$. The quadrupole $({\ell },| m| )=(2,2)$ mode remains dominant during most of the inspiral, while other, higher-order multipoles become comparable near the merger and ringdown stages (see, e.g., Pan et al. 2011; Calderón Bustillo et al. 2016). Higher-order multipoles have been found to be important for high-mass BBH mergers in advanced detector observations (Pan et al. 2011; Varma et al. 2014; Graff et al. 2015; Varma & Ajith 2017; Calderón Bustillo et al. 2017, 2018; Mehta et al. 2017; Chatziioannou et al. 2019; Kumar Mehta et al. 2019); their presence was also crucial in reducing uncertainties in the mass ratio and component spins of the recently reported unequal-mass binary mergers GW190412 (Abbott et al. 2020a) and GW190814 (Abbott et al. 2020e). To quantify the impact and evidence of higher-order multipoles on GW190521, we computed the posterior distribution of our preferred model, NRSur PHM (Varma et al. 2019), using the RIFT pipeline (Lange et al. 2018) for three different combinations of modes: ℓ = 2, ℓ ≤ 3, and ℓ ≤ 4.

GR predicts radiation from BBH mergers at all multipoles, with amplitudes strongly dependent on the binary mass ratio and orbital inclination angle (Graff et al. 2015; London et al. 2018); evidence of their presence consistent with GR is presented in Abbott et al. (2020a).

We find that the omission of higher-order multipoles leads to broader posterior distributions for some parameters of GW190521, especially the binary orbital inclination angle θ_JN. The higher-order multipoles enable better constraints of the binary inclination angle, which is coupled to the source-frame mass estimates through the luminosity distance D_L, and therefore the redshift. Inclusion of these higher-order multipoles results in significantly larger values for these parameters. As shown in Figure 7, this change in distance directly impacts our inference of the primary BH mass of GW190521. For our preferred model, if we only included ℓ = 2, the posterior distribution of primary BH mass extends to m₁ = 137 M_⊙ at the 90% credible interval. As we include higher-order multipoles up to ℓ = 4, the masses are better constrained and we find that m₁ cannot be greater than 103 M_⊙ at the 90% credible interval. We find that the higher-order multipoles have marginal impact on constraining the effective precession spin parameter χ_p.

Zoom In Zoom Out Reset image size

Although the released version of NRSur PHM (Varma et al. 2019) does not include multipoles with ℓ > 4, we have extended it to include all multipoles with ℓ = 5 and find that the parameter estimation results are nearly identical to those for ℓ ≤ 4; the ℓ = 5 multipoles are found to be quantitatively negligible for GW190521.

While higher-order multipoles influence our inference of the source properties, the observation of GW190521 does not offer evidence for such multipoles. The log₁₀-Bayes factor for our preferred model with and without higher-order multipoles is −0.38, implying that the data marginally disfavor their presence, but there is no statistically significant evidence for their absence. As noted above, this gives important information about the orbital inclination of the binary, as higher-order multipoles are suppressed for face-on binary orbits. Since higher-order multipoles are a firm prediction of GR and significantly affect parameter estimation, they must be included when modeling GW190521.

We compared the data from LIGO and Virgo at the time of GW190521 directly with 3459 NR simulations of BBH coalescence (Jani et al. 2016; Boyle et al. 2019; Healy et al. 2019). These NR simulations are the most accurate representation available of the strong-field dynamics near merger and the radiated higher-order multipoles. As was done for previous direct comparisons of NR simulations with LIGO–Virgo-detected events (Abbott et al. 2016a, 2019i; Healy et al. 2018), the likelihood of the data was calculated for each waveform derived from NR simulations. Using the RIFT pipeline (Lange et al. 2018) to evaluate these comparisons and interpolate them over all intrinsic parameters, we deduce a posterior distribution for all detector-frame quantities: redshifted mass, mass ratio, and both component spins.

As shown in Figure 8, we find that the best-matching (highest-likelihood) NR waveforms and the derived posterior distribution are consistent with the mass ratio of GW190521 being near unity. We also find that the analysis favors NR waveforms with χ_eff ∼ 0 and χ_p ∼ 0.6, further suggesting that GW190521 is a precessing binary. The agreement between our preferred model NRSur PHM and the NR simulations provides an independent check on the inferred source properties of GW190521.

Zoom In Zoom Out Reset image size

We first test the consistency of the templates used for parameter estimation with the observed signal by subtracting the maximum likelihood NRSur PHM waveform from the data and studying the residual (Abbott et al. 2016c, 2020f). The residual data are analyzed by the BayesWave algorithm, which simultaneously fits models for the noise PSD and short (<1 s) transient signals coherent across the detector network (Cornish & Littenberg 2015). BayesWave does not assume a particular waveform morphology, instead using a linear combination of wavelets to fit excess coherent features in the data of arbitrary shape.

The metric used for quantifying the coherent residual, S/N₉₀, is the upper 90% credible bound on the posterior distribution function of the recovered S/N of the wavelet reconstruction (Abbott et al. 2016c). The significance of the recovered S/N₉₀ value is empirically measured by repeating the analysis on 195 randomly chosen off-source times drawn from 4096 s of detector data surrounding the event time, which serves as the background estimate for S/N₉₀. For GW190521, we find S/N₉₀ ∼ 6.34, which is consistent with expectations for typical LIGO–Virgo noise, resulting in a p-value of ∼0.26 when compared to the times analyzed immediately surrounding the merger time. Since the residual is consistent with noise, we find that the best-fit waveform interpolated from numerical solutions of GR, as used in our parameter estimation analyses, is consistent with the observed signal within the measurement capability of the detectors.

In this section, we describe a number of studies that explore the nature of the remnant compact object and also evaluate the consistency of the observed signal with waveform models for a quasi-circular BBH merger in GR. These are extensions to similar tests applied to previous LIGO–Virgo BBH events (Abbott et al. 2016c, 2019g). As in the case for those previous tests, the high mass of GW190521 makes it possible to study the quasi-normal oscillations of the remnant BH approaching a stationary state (ringdown), as encoded in the last few cycles of the GW signal (see Vishveshwara 1970; Buonanno et al. 2007; Berti et al. 2009, 2018, and references therein). Parameter estimates obtained from the ringdown studies presented in this section are solely extracted from the properties of the remnant. For these studies, we consider both a model that makes no assumptions on the process leading to the signal formation and waveform templates specifically modeling the remnant of a quasi-circular BBH merger. The results are thus robust against systematic uncertainty due to the possible neglect of eccentricity or other physical effects in modeling the system; for further discussion of such effects see Section 6.1.

We model the ringdown as a set of damped sinusoids and measure the properties of the signal using pyRing, a time domain analysis framework (Carullo et al. 2019; Isi et al. 2019). Within the analysis, the beginning of the ringdown-dominated portion of the signal is marked by a fixed time t₀ (see Supplement of Isi et al. 2019), reported with respect to an estimate of the peak time of the complex strain at each detector. The peak time is taken from the posterior median inferred from the NRSur PHM model and yields a GPS time for LIGO Hanford ${t}_{\mathrm{peak}}^{{\rm{H}}}={1242442967.4306}_{-0.0106}^{+0.0067}\,{\rm{s}};$ times in other detectors are computed by assuming a fixed sky position. The duration of the analysis segment after the start time in each detector is 0.1 s. We consider three distinct models of the ringdown: for each, we use the CPNest Bayesian nested sampling algorithm (Del Pozzo & Veitch 2020) to compute posterior distributions on the model parameters.

The first model consists of a single damped sinusoid of arbitrary complex frequency, amplitude, and phase, with no assumptions on the nature of the remnant object. The left panel of Figure 9 shows the resulting 90% credible regions for the posterior probability distribution of the redshifted, detector-frame frequency and damping time, assuming uniform priors on these two parameters (solid curves). Since there is uncertainty in the time at which the single damped sinusoid model may become valid, we perform this analysis by truncating the data at different start times, as in Kamaretsos et al. (2012) and Abbott et al. (2016c): we choose start times ${\rm{\Delta }}{t}_{0}\equiv {t}_{0}-{t}_{\mathrm{peak}}^{{\rm{H}}}\simeq 6.4$, 12.7, and 19.1 ms after the reference time. In units of the redshifted remnant mass ${M}_{{\rm{f}}}^{\det }=(1+z){M}_{{\rm{f}}}=258.3$ M_⊙ (median estimate from NRSur PHM) these times correspond to ${\rm{\Delta }}{t}_{0}=5,10,15\,{{GM}}_{{\rm{f}}}^{\det }/{c}^{3}$. Truncating the data at later times yields uninformative posteriors due to the decrease in S/N. Values of ${\rm{\Delta }}{t}_{0}\geqslant 10\,{{GM}}_{{\rm{f}}}^{\det }/{c}^{3}$ are consistent with the range where numerical simulations predict that the fundamental mode should be dominant for this source. The systematic uncertainty due to the choice of the NRSur PHM model in determining t_peak as opposed to other waveform models considered in this paper is negligible with respect to the statistical uncertainty.

Zoom In Zoom Out Reset image size

At a start time of Δt₀ = 12.7 ms (corresponding to $10\,{{GM}}_{{\rm{f}}}^{\det }/{c}^{3}$), we measure a ringdown frequency f = ${66.0}_{-3.0}^{+4.0}$ Hz and damping time τ = ${19.0}_{-7.0}^{+9.0}$ ms. As the start time increases from 6.4 to 19.1 ms, we observe convergence toward the predicted value of the least damped mode frequency, together with a broadening of the posterior due to the decreasing S/N of the ringdown. Similar behavior is observed for the other models we will present in the remainder of the section.

The second model is built from the superposition of a set of damped sinusoids with arbitrary amplitudes and phases, but with complex frequencies determined by the remnant mass M_f and dimensionless spin χ_f, as predicted by perturbation theory (Berti et al. 2006). We include up to the n = 0, 1, 2 overtones of the ℓ = m = 2 ringdown mode of a Kerr BH starting at the peak of the complex strain (Buonanno et al. 2007; Giesler et al. 2019; Ota & Chirenti 2020); see also Bhagwat et al. (2020) and Okounkova (2020) for discussions concerning the interpretation of a linearized approximation starting at the peak of GW emission. The redshifted, detector-frame remnant mass and spin obtained from this waveform model are shown in the right panel of Figure 9, assuming a uniform prior over these two parameters. We find ${M}_{{\rm{f}}}^{\det }={252.0}_{-64.0}^{+63.0}$ M_⊙, χ_f = ${0.65}_{-0.48}^{+0.22}$ taking Δt₀ = 12.7 ms and including only the fundamental Kerr (ℓ = 2, m = 2, n = 0) mode, and ${M}_{{\rm{f}}}^{\det }\,={276.0}_{-50.0}^{+51.0}$ M_⊙, χ_f = ${0.74}_{-0.3}^{+0.15}$ taking Δt₀ = 0 ms and including overtones up to n = 2.

Finally, the third model (Kerr HMs) consists of a set of damped sinusoids corresponding to all fundamental modes (i.e., without inclusion of overtones) of a Kerr BH up to ℓ = 4 and m = ℓ, ℓ − 1, including spherical−spheroidal harmonic mixing (London 2018). Complex frequencies are predicted as a function of the remnant mass M_f and dimensionless spin χ_f, while amplitudes and phases are calibrated on NR simulations of nonprecessing BBH mergers. With this model we measure ${M}_{{\rm{f}}}^{\det }={296.0}_{-58.0}^{+52.0}$ M_⊙, χ_f = ${0.79}_{-0.3}^{+0.12}$ taking Δt₀ = 12.7 ms; the full probability distribution is shown in the right panel of Figure 9.

In Figure 9, we compare the ringdown measurements to the posterior credible regions for the remnant parameters obtained through NR-calibrated fits from the initial binary parameters, as described in Section 2.1, from the three different full waveform models discussed above. The posterior of the ringdown analyses is consistent at the 90% credible level with the full-signal analyses. Furthermore, despite the different physical content, both models that include higher multipoles or overtones obtain measurements of remnant parameters consistent with the single-mode analysis estimates; a Bayes factor computation also does not find strong evidence in favor of the presence of higher multipoles or overtones.

A parameterized test of gravitational waveform generation (Blanchet & Sathyaprakash 1995; Mishra et al. 2010; Li et al. 2012a, 2012b; Agathos et al. 2014) based on the Phenom PHM waveform model also does not reveal inconsistencies with GR predictions. Full details will be provided in an upcoming paper.

The physical processes leading to PI and PPI are well known, and a robust formalism to model them has been developed (e.g., Fowler & Hoyle 1964; Barkat et al. 1967; Ober et al. 1983; Bond et al. 1984; Heger et al. 2003; Woosley et al. 2007; Chen et al. 2014; Yoshida et al. 2016; Woosley 2017, 2019; Marchant et al. 2019). On the other hand, there are still uncertainties on the minimum helium-core mass for a star to undergo PPI: Woosley (2017) indicate M_He,min ≈ 32 M_⊙, while Leung et al. (2019) suggest that this value should be raised to M_He,min ≈ 40 M_⊙. This difference might translate into a significant difference in the maximum BH mass.

In addition to this, there are critical uncertainties on other physical ingredients that combine with PI and PPI to shape the mass spectrum of BHs (e.g., Belczynski et al. 2016, 2020; Stevenson et al. 2019). Recently, Farmer et al. (2019) and Renzo et al. (2020) have investigated the main sources of uncertainty on the lower edge of the PI mass gap, by modeling naked helium cores. The impact of time-dependent convection on the location of the lower edge of the PI mass gap is found to be Δm ≈ 5 M_⊙ (Renzo et al. 2020). Variations in the treatment of convective mixing and neutrino physics are found to have small effects on the maximum BH mass (with a mass variation Δm ≈ 1–2 M_⊙), while different stellar metallicity and wind mass-loss prescriptions have more conspicuous repercussions (Δm ≈ 4 M_⊙). Most importantly, the uncertainties on nuclear reaction rates have a dramatic impact on the maximum BH mass, which can change by Δm ≈ 16 M_⊙, corresponding to 30%–40% of the maximum BH mass. Most of the differences come from the ¹²C(α, γ)¹⁶O reaction (Farmer et al. 2019; see also Takahashi et al. 2018). This result is particularly important if we consider that Farmer et al. (2019) explore only 1σ uncertainties (see, e.g., deBoer et al. 2017 for a recent review) and vary only a few nuclear reaction rates.

The maximum BH mass estimated by Farmer et al. (2019) is ≈56 M_⊙; however, they model pure helium cores, without hydrogen envelopes. Mapelli et al. (2020) show that accounting for the collapse of a residual hydrogen envelope can increase the maximum BH mass by ≈20 M_⊙ in the case of metal-poor (Z ≤ 0.0003), slowly rotating progenitor stars. Hence, the final fate of the hydrogen envelope is an additional source of uncertainty, with an impact of Δm ≈ 20 M_⊙ on the maximum BH mass.

While mass transfer in close binaries is expected to remove most of the hydrogen envelope, metal-poor (Z ≤ 3 × 10⁻⁴) single stars with zero-age main-sequence (ZAMS) mass m_ZAMS ≲ 70 M_⊙ are expected to retain a significant fraction of their hydrogen envelope. If most of this hydrogen envelope collapses into the final compact object, the maximum BH mass is ∼60–65 M_⊙, while a naked helium core can produce BHs up to ∼45 M_⊙ (assuming median values for the nuclear reaction rates; Mapelli et al. 2020).

In summary, many sources of uncertainty affect our knowledge of the PI mass gap, the two principal ones being nuclear reaction rates (Δm ≈ 16 M_⊙) and the collapse of the hydrogen envelope (Δm ≈ 20 M_⊙). The combination of the main uncertainties has yet to be studied; therefore, it is not clear whether their effects sum up, as needed to explain the primary mass of GW190521. Based on these considerations, ∼60–65 M_⊙ is probably a conservative lower limit to the edge of the mass gap.

In this section, we discuss general properties of hierarchical mergers, i.e., mergers involving one or two second-generation BHs; we then present a Bayesian estimate of the relative probabilities (odds) that GW190521 is the product of a hierarchical BH merger, as opposed to a merger of first-generation (hereafter 1g) BHs. The masses of 2g BHs can fall in the PI gap: the merger remnant of GW170729 (Abbott et al. 2019i) is a previously reported example. A 2g BH is alone at birth, unless it was a member of a triple system that remains bound after the formation of the BH. However, a single 2g BH might acquire a companion through dynamical exchanges, if it forms in a dense cluster (Miller & Hamilton 2002; Colpi et al. 2003; Rodriguez et al. 2019), or through migration-mediated interactions in an AGN disk (McKernan et al. 2012; Bartos et al. 2017).

When formed, a 2g BH is subject to a relativistic gravitational recoil velocity (kick), which can eject it from its birthplace (Merritt et al. 2004); see Figure 4. For example, if a 2g BH forms in a dense star cluster, it might be retained in the cluster and subsequently acquire a new companion only if the local escape velocity is larger than ≳50 km s⁻¹ (Holley-Bockelmann et al. 2008; Antonini & Rasio 2016; O'Leary et al. 2016; Gerosa & Berti 2019). In this scenario, the companion might be a 1g BH or, less likely, a 2g BH (Rodriguez et al. 2019).

The relativistic kick at birth strongly depends on BH spin: maximally spinning BHs with spins in the orbital plane and counteraligned receive the largest kicks (e.g., Campanelli et al. 2007a, 2007b). Under the simplified assumption that all 1g BHs are born with spin χ = 0, then ∼60% of 2g BHs are retained in a typical globular cluster, while if all 1g BH have χ = 0.5, then less than 3% of merger products are retained (Rodriguez et al. 2019). For these reasons, hierarchical mergers are expected to be more likely in nuclear star clusters (e.g., Antonini & Rasio 2016; Arca Sedda et al. 2020), which have the highest escape velocity, rather than globular clusters or smaller stellar systems, such as open clusters and young star clusters.

Binaries with 2g BHs have several distinctive properties. The mass of the primary can be significantly larger than the lower boundary of the PI mass gap, although the exact distribution strongly depends on the 1g mass function (Chatziioannou et al. 2019; Gerosa & Berti 2019; Rodriguez et al. 2019; Doctor et al. 2020). Concerning the predicted mass ratios, 1g+2g mergers tend to have more unequal mass ratios than 2g+2g mergers (Rodriguez et al. 2019). The preference of GW190521 for a mass ratio q ∼ 1 suggests that this system is more likely consistent with a 2g+2g merger than with a 1g+2g merger. Indeed, it is also conceivable for the system to be a merger of 1g BHs where one component was above the PI mass gap, though this is disfavored by the preference for near-equal masses (see Figure 8).

Spin measurements are a distinguishing feature of 2g BHs, because merger remnants are, on average, rapidly rotating with χ ∼ 0.7 (e.g., Baker et al. 2004; Buonanno et al. 2008; Lousto et al. 2010; Hofmann et al. 2016; Fishbach et al. 2017). Since they acquire companions through dynamical exchanges, spins are expected to be isotropically distributed with respect to the BBH orbital angular momentum (e.g., Mandel & O'Shaughnessy 2010; Rodriguez et al. 2016). A direct consequence of spatial isotropy is that the predicted distribution of effective spins is broad and symmetric about χ_eff = 0 (Farr et al. 2018; Gerosa & Berti 2019; Rodriguez et al. 2019; Kimball et al. 2020a). This prediction is broadly consistent with the range of χ_eff and relatively large χ_p values inferred for GW190521, reported in Section 2.2.

Finally, to estimate the predicted merger rate of PI BBHs from the hierarchical scenario, we should take into account that these mergers come from all kinds of star clusters with v_esc > 50 km s⁻¹, hence both massive globular clusters and nuclear star clusters (or AGN disks). Considering only globular clusters and assuming that 1g BHs have zero spin, Rodriguez et al. (2019) find that ∼3% of all mergers from globular clusters at redshift z < 1 have a component greater than 55 M_⊙, corresponding to 7% of the detectable BBHs from globular clusters.

5.2.1. Bayesian Analysis of Hierarchical Formation in Globular Clusters

We apply the method of Kimball et al. (2020b) to recover posteriors over the relative rates of 1g+2g and 2g+2g mergers to 1g+1g mergers, as well as the odds ratios that the GW190521 binary is of 1g+2g or 2g+2g origin, using our estimates of source component masses and spins within a physical model of the hierarchical merger process adapted to globular clusters. This analysis illustrates the many physical parameters and associated uncertainties that enter any inference on the origin of massive BBH.

We perform hierarchical Bayesian inference on a population model containing 1g+1g, 1g+2g, and 2g+2g mergers. For the first-generation BBH population, we adopt Model C in Abbott et al. (2019h): the mass distribution is a mixture of a truncated power law and a high-mass Gaussian component (Talbot & Thrane 2018), while the component spins follow nonsingular beta distributions with isotropic orientations. We also apply a theoretically motivated prior to the upper cutoff of the power-law mass component m_max to reflect expectations from stellar collapse dynamics, in particular PI (see Section 5.1). Considering uncertainties due to nuclear reaction rates (Farmer et al. 2019) and fallback dynamics (Mapelli et al. 2020), we take a prior Gaussian distribution over m_max with mean 50 M_⊙ and standard deviation 10 M_⊙. As mentioned above, we might also consider a 1g component above the PI mass gap, i.e., of ∼120 M_⊙; however, this would add complexity to the analysis while being unlikely to significantly alter the main outcome.

For simplicity, all 1g+1g mergers are here assumed to occur in environments that could potentially lead to subsequent mergers (Doctor et al. 2020). The 1g+2g and 2g+2g populations are formed by applying transfer functions from Kimball et al. (2020b) to the 1g+1g distributions, motivated by simulations reported in Rodriguez et al. (2019). The populations are combined into a mixture model by calculating branching ratios for the fraction of 1g+1g merger remnants that are retained in their cluster environments after receiving relativistic kicks at merger (calculated with the precession package of Gerosa & Kesden 2016; Gerosa 2016). The retention probability depends on the cluster escape velocity: as in Kimball et al. (2020b), we take a nominal cluster mass of 5 × 10⁵ M_⊙ with a Plummer radius of 1 pc (Harris 1996; Kremer et al. 2019); to quantify dependence on the escape velocity, we also consider a higher mass cluster of 10⁸ M_⊙ representative of nuclear cluster environments.

The overall retention fraction is also highly sensitive to the 1g+1g spin magnitude distribution. Fuller & Ma (2019) find that a significant fraction of stellar BHs should form with near-zero natal spins; for BBHs formed in clusters via dynamical interaction we may also neglect any possible effect of strong stellar binary interactions on 1g BH spins. Thus, we also perform an analysis with an expanded model that includes an additional subpopulation of 1g BHs with zero spin, making up a fraction λ₀ of 1g+1g binary components. We quote results both with and without this zero-spin channel.

We apply our mixture models to a set of BBH observations consisting of GW190521 and all BBH detections from O1 and O2, performing inference on model hyperparameters using GWPopulation (Talbot et al. 2019) and Bilby (Ashton et al. 2019). In Figure 11, we plot the resulting marginalized Bayes factor of a given BBH being of xg+2g versus 1g+1g origin as a function of primary component mass and spin, together with the 90% and 68% credible regions for GW190521. Overall Bayes factors including dependence on both component masses and spins generally favor 2g+2g origin; the parameter estimation likelihood strongly favors near-equal component masses; thus, 1g+2g origin implying significantly unequal masses is disfavored by the data.

Zoom In Zoom Out Reset image size

Figure 12 (orange contours) shows the relative rates of 1g+2g and 2g+2g mergers for the model without a 1g zero-spin subpopulation. We find ${\mathrm{log}}_{10}\left({{ \mathcal R }}_{1{\rm{g}}+2{\rm{g}}}/{{ \mathcal R }}_{1{\rm{g}}+1{\rm{g}}}\right)$ and ${\mathrm{log}}_{10}\left({{ \mathcal R }}_{2{\rm{g}}+2{\rm{g}}}/{{ \mathcal R }}_{1{\rm{g}}+1{\rm{g}}}\right)$ to be $-{2.6}_{-1.1}^{+0.9}$ and $-{5.5}_{-2.2}^{+1.8}$, respectively. Applying the priors given by the relative rates and marginalizing over system parameters and population hyperparameters, we calculate the odds ratio that the source of GW190521 is a hierarchical merger, versus a 1g+1g merger. Our results, summarized in the top half of Table 2, are dependent on the waveform chosen to produce the parameter estimation samples, as well as assumptions on the 1g+1g spin and mass ratio distributions. Under these assumptions, we find that GW190521 is most likely to be of 1g+1g origin, because the merger remnants of BBHs with large component spins are often subject to kicks that eject them from low-mass clusters. Accounting for a possible 1g BH component above the mass gap, implying a significantly unequal-mass merger, would be unlikely to alter this conclusion.

Zoom In Zoom Out Reset image size

Table 2.  Odds Ratios for the Source of GW190521 Being a xg+2g Merger versus 1g+1g Merger, in the Globular Cluster Models of Section 5.2.1

No Zero-spin Channel	NRSur PHM	Phenom PHM
P(1g+2g):P(1g+1g)	0.28:1	0.28:1
P(2g+2g):P(1g+1g)	0.07:1	0.03:1
With Zero-spin Channel	NRSur PHM	Phenom PHM
P(1g+2g):P(1g+1g)	0.86:1	0.81:1
P(2g+2g):P(1g+1g)	0.68:1	0.23:1

Download table as:  ASCIITypeset image

Applying the model with a 1g zero-spin subpopulation yields higher retention fractions, due to lower merger recoil velocities, and higher relative hierarchical merger rates, shown in Figure 12 (blue contours): ${\mathrm{log}}_{10}\left({{ \mathcal R }}_{1g+2g}/{{ \mathcal R }}_{1{\rm{g}}+1{\rm{g}}}\right)$ and ${\mathrm{log}}_{10}\left({{ \mathcal R }}_{2{\rm{g}}+2{\rm{g}}}/{{ \mathcal R }}_{1{\rm{g}}+1{\rm{g}}}\right)$ become $-{2.0}_{-1.2}^{+1.0}$ and $-{4.4}_{-2.4}^{+2.0}$, respectively. When we include the zero-spin formation channel, the odds ratios for hierarchical origin, summarized in the bottom half of Table 2, increase by factors of ∼2–10, though we still favor GW190521 being due to a 1g+1g merger, for both the NRSur PHM and Phenom PHM parameter estimates. In this case, the estimated fraction λ₀ of merging 1g BHs from the zero-spin channel is in the range of ∼0.01–0.47.

We also consider a higher cluster mass of 10⁸ M_⊙, which may be representative of nuclear cluster environments: the higher cluster mass increases the odds of hierarchical origin by 3–4 orders of magnitude, for both models with and without a zero-spin formation channel. For AGN environments essentially all 1g merger products are retained; however, quantitative assessment would require more detailed modeling of BH mergers in AGNs beyond their increased escape velocity, and the 1g+2g and 2g+2g population models we present here are tuned to globular cluster models.

In summary, the probability that GW190521 is due to a hierarchical merger in a stellar cluster is strongly dependent on the properties of 1g BHs in such environments, primarily on their mass and spin distributions and on the cluster escape velocity. With significantly larger event samples it may be possible to disentangle the different model parameters; thus, similar, high-mass BBH mergers constitute a future laboratory for BH populations and dynamics in cluster environments.

If a star grows an oversized hydrogen envelope with respect to its helium core, it might directly collapse to a BH with mass ∼60–100 M_⊙ without entering the PI/PPI regime (Bouffanais et al. 2019; Di Carlo et al. 2019; Spera et al. 2019). This scenario assumes that most of the hydrogen envelope collapses to a BH (see Sukhbold et al. 2016 for a discussion of the uncertainties). For a star to develop this oversized hydrogen envelope, one or more mergers between a helium-core giant and a main-sequence companion are required (Di Carlo et al. 2020). As in the hierarchical merger scenario, a BH born from a stellar merger is alone at birth, unless it was a member of a triple system that remains bound. In the field the BH remains alone, while in a young star cluster, a globular cluster, or a nuclear star cluster the BH can acquire a companion through exchanges. A further variation of this model is that GW190521 formed in a triple or multiple stellar system. The inner binary star of the triple might have merged, producing the primary BH via the stellar merger mechanism described above. If the outer member is massive enough, it might then collapse, becoming the secondary BH. The most critical aspect of a triple origin for GW190521 is the orbit of the tertiary body.

The primary mass distribution predicted by the stellar merger model scales approximately as m⁻⁵ with m > 60 M_⊙. The expected mass ratio is most likely q ∼ 0.4–0.6, but all mass ratios between q ∼ 0.04 and q ∼ 1 are possible (Di Carlo et al. 2020). Hence, the primary mass and the mass ratio predicted by this scenario are compatible with the properties of GW190521.

If the angular momentum of the core of the massive progenitor is not dissipated efficiently, this scenario predicts a high primary spin. Since this is a dynamical formation scenario, spin orientations are expected to be isotropically distributed. As we already discussed for the hierarchical scenario (Section 5.2), an isotropic spin distribution is consistent with the fact that GW190521 has a low χ_eff and a relatively large χ_p.

Finally, Di Carlo et al. (2019) predict that up to ∼2% of all BBH mergers from young star clusters involve BHs in the PI mass gap born from (multiple) stellar mergers. These represent ≲10% of all detectable mergers from star clusters. As noted in Jani & Loeb (2020), ≲0.8% of all massive stars contribute to a BBH merger population in the PI mass gap.

In summary, the primary mass, mass ratio, effective spin, and precession spin parameters of GW190521 are consistent with the stellar merger scenario in star clusters. The key difference between the hierarchical merger scenario and the stellar merger scenario is that the latter does not imply relativistic kicks at birth; hence, it might be at least one order of magnitude more common in star clusters. On the other hand, the details of the stellar merger scenario depend on delicate assumptions about stellar mergers (e.g., that most mass is retained by the merger product) and massive star evolution (e.g., stellar rotation).

The nucleus of an active galaxy might harbor tens of thousands of stellar-mass BHs that moved into the innermost parsec owing to mass segregation (Bahcall & Wolf 1977; Morris 1993; Miralda-Escudé & Gould 2000; Antonini 2014; Hailey et al. 2018; Generozov et al. 2018). In this dense gaseous environment, BH orbits are efficiently torqued by gas drag until they align with the AGN disk (Bartos et al. 2017; McKernan et al. 2018). Once in the disk, BHs can accrete gas, and they are expected to acquire companions and efficiently merge with them as an effect of gas torques (McKernan et al. 2012, 2014, 2020; Bellovary et al. 2016; Bartos et al. 2017). When these BHs merge with other BHs, there is a high chance that they are not ejected, because of the high escape velocity in galactic nuclei (thousands of kilometers per second, given the proximity to an SMBH). Hence, AGN disks might easily harbor 2g BHs (Yang et al. 2019a).

As a consequence of the joint contribution of gas accretion, BBH mergers, and high escape velocities, a fraction of BHs in AGN disks are expected to have masses in the PI mass gap or even in the IMBH regime (McKernan et al. 2012). Recent work suggests that they have a mass function similar to a power law, but significantly flatter than field BHs (Yang et al. 2019b), or even reminiscent of a broken power law (McKernan et al. 2018; Yang et al. 2019a). The possible mass ratios are highly uncertain, though most mergers have mass ratios less extreme than 10:1 (McKernan et al. 2018; Yang et al. 2019a).

As to the spin magnitudes, we expect large spins (χ ∼ 0.7) if BHs in AGN disks are 2g BHs. Since BH orbits tend to align with the AGN disk, there might be some preferential alignment in the spins of BBHs (McKernan et al. 2012). Unlike the other possible scenarios, BBH mergers in an AGN disk may have an associated electromagnetic counterpart (an ultraviolet flare; McKernan et al. 2019). No prompt electromagnetic counterpart has been reported for GW190521 (e.g., Barthelmy et al. 2019; Casentini et al. 2019; Doyle et al. 2019a, 2019b; Lipunov et al. 2019; Negoro et al. 2019; Turpin et al. 2019a, 2019b; Watson et al. 2019), but see Graham et al. (2020), which appeared during the revision of this work, for a possible candidate optical counterpart.

Finally, the fraction of BBHs in the PI mass gap born in AGN disks might be large: Yang et al. (2019a) suggest that 40% of all AGN-assisted mergers detected by LIGO–Virgo will include a BH with mass ≳50 M_⊙; on the other hand, the overall contribution of AGN disks to the BBH merger rate might be low. Most works suggest a total merger rate of ∼1–10 Gpc⁻³ yr⁻¹ (e.g., Bartos et al. 2017; Stone et al. 2017; Yang et al. 2019b; Arca Sedda 2020), while McKernan et al. (2018) and Tagawa et al. (2020) try to quantify additional uncertainties and end up with merger rates of ∼10⁻³–10⁴ Gpc⁻³ yr⁻¹ and ∼0.02–60 Gpc⁻³ yr⁻¹, respectively.

The waveform observed in the data may also be consistent with the merger of a BBH with nonzero orbital eccentricity. The short duration of the signal makes it difficult to distinguish the amplitude modulation associated with precession from that due to eccentric orbits, or even head-on collisions (mergers that happen immediately at closest approach; Calderón Bustillo et al. 2020). More quantitative evaluation of this issue requires further development of accurate and computationally efficient eccentric waveforms and is thus deferred to future work. Efforts to efficiently detect and estimate the parameters for eccentric compact binaries are under development (Lower et al. 2018; Romero-Shaw et al. 2019; Lenon et al. 2020).

From an astrophysical perspective, BBHs with eccentricity e ≳ 0.1 at merger are deemed to be almost impossible from isolated binary evolution but might account for ∼1% of all BBH mergers in globular clusters and other dense star clusters (Gültekin et al. 2006; Samsing et al. 2014, 2018; Samsing & Ramirez-Ruiz 2017; Rodriguez et al. 2018; Samsing 2018; Samsing & D'Orazio 2018; Zevin et al. 2019; Fragione et al. 2020). Head-on collisions are even more rare owing to the small geometric cross section of individual BHs (Samsing et al. 2014).

The high concentration of BHs in nuclear star clusters might produce a population of extremely eccentric BBHs from single–single BH capture by GW radiation (O'Leary et al. 2009; Gondán et al. 2018), although the cross section for this process is orders of magnitude lower than the cross section for a binary—single encounter (Samsing et al. 2014). Finally, Kozai–Lidov resonances (Kozai 1962; Lidov 1962) in triple systems might also trigger eccentric BBHs in the field (Antognini et al. 2014; Antonini et al. 2017; Silsbee & Tremaine 2017) or in dense star clusters (Wen 2003; Antonini & Perets 2012; Antonini et al. 2016; Kimpson et al. 2016; Hoang et al. 2018). Considering that ∼25% of massive stars are members of triple systems (Sana et al. 2014), Antonini et al. (2017) estimate a merger rate of ∼0.3–2.5 Gpc⁻³ yr⁻¹ from isolated resonating triple systems, up to ∼5% of which retain high eccentricity at merger.

Most of the formation channels that can explain the mass of the primary BH in GW190521 are based on dynamical encounters in dense star clusters and allow for the formation of eccentric BBHs. Hence, even if eccentric BBH mergers are estimated to be extremely rare, we cannot exclude that the source binary of GW190521 has nonzero eccentricity.

Another possible explanation for the apparent inconsistency of GW190521 with formation from stellar collapse is gravitational lensing of the signal by galaxies or galaxy clusters (Broadhurst et al. 2018; Li et al. 2018; Ng et al. 2018; Oguri 2018; Smith et al. 2018). If GW190521 is a strongly lensed signal, it will receive a magnification μ defined such that GW amplitude is increased by a factor μ^1/2 relative to the unlensed case. For a given observed GW signal, the luminosity distance to the source, and thus its redshift, may then be much larger than in the absence of lensing. Strong lensing is expected to be relatively rare at current detector sensitivities, with 1 in every 100–1000 detected events strongly lensed by individual galaxies (Holz & Wald 1998; Li et al. 2018; Ng et al. 2018; Oguri 2018) and a similarly low rate for lensing by galaxy clusters (Smith et al. 2018; Oguri 2019).

Under the strong-lensing hypothesis ${{ \mathcal H }}_{L}$, we assume that the event comes from a BBH population of stellar collapse origin and infer its magnification, redshift, and component masses. The posterior distribution of the parameters is (Pang et al. 2020)

$$\begin{eqnarray}&&p(\mu ,\vartheta | d,{{ \mathcal H }}_{L})\propto p(d| \vartheta )p(\vartheta | \mu ,{{ \mathcal H }}_{L})p(\mu | {{ \mathcal H }}_{L}),\end{eqnarray} \tag{ 1 }$$

where ϑ are the apparent parameters of the waveform received at the detector, which differ from the source-frame parameters owing to effects of redshift and lensing. For the prior over component masses and redshift of the source $p({m}_{1},{m}_{2},z| {{ \mathcal H }}_{L})$, we take a binary BH population model from those used in Abbott et al. (2019h) for O1 and O2 observations, with fixed population parameters λ = 0, α = 1, β_q = 0, m_min = 5 M_⊙, multiplied by the optical depth of lensing by galaxies τ(z) = 4.17 × 10⁻⁶ (D_c(z)/Gpc)³ (Haris et al. 2018), where D_c is the comoving distance. For the primary mass upper cutoff m_max we consider two different values to account for uncertainties in the edge of the PISN gap, m_max = (50, 65) M_⊙. We use the lensing prior $p(\mu | {{ \mathcal H }}_{L})\propto {\mu }^{-3}$ (Blandford & Narayan 1986) with a lower limit μ > 2 appropriate to strong lensing. We adopt the NRSur PHM waveform model. The resulting magnification estimate is mostly sensitive to the measurement of the component masses.

We find that the required magnification under the lensing hypothesis, taking m_max = 50 M_⊙ (m_max = 65 M_⊙), is ${13}_{-8}^{+54}$ (${6}_{-4}^{+28}$), with source-frame masses ${43}_{-16}^{+6}$ M_⊙ (${55}_{-22}^{+9}$ M_⊙) and ${35}_{-15}^{+11}$ M_⊙ (${44}_{-19}^{+14}$ M_⊙) at redshift ${2.5}_{-0.7}^{+2.1}$ (${1.8}_{-0.6}^{+1.7}$). At these redshifts, the lensing optical depth is low: ${9}_{-3}^{+10}\times {10}^{-4}$ (${5}_{-3}^{+9}\,\times {10}^{-4}$).

One possible signature of strong lensing would be multiple GW images, which may give rise to two events occurring closer in time than expected for a Poisson process of BBH mergers, and with consistent sky localization and intrinsic parameters. Another candidate GW signal, S190521r (Abbott et al. 2019c), was reported 4.6 hr after GW190521; however, the two events' sky localizations strongly disfavor lensing, showing no overlap (Haris et al. 2018; Hannuksela et al. 2019; Singer et al. 2019). Any lensed counterpart image could, though, have been too weak to be confidently detected and may have arrived at a larger time separation or when the detectors were not operating. Moreover, in the case of galaxy cluster lensing, the counterpart could appear at a separation of years in time (Smith et al. 2018).

Given the low expected lensing rate and optical depth and the absence of an identifiable multi-image counterpart close to GW190521, our current analyses find no evidence in favor of the strong-lensing hypothesis. Future analyses using subthreshold searches and multi-image searches on all event pairs may yield better constraints on strong lensing (Haris et al. 2018; Hannuksela et al. 2019; Li et al. 2019; McIsaac et al. 2019).

Primordial BHs (PBHs; Carr & Hawking 1974; Khlopov 2010) are thought to be formed from collapse of dark matter overdensities in the very early universe (at redshifts z > 20, i.e., before the formation of the first stars) and may account for a nontrivial fraction of the density of the universe (Carr et al. 2016; Clesse & García-Bellido 2017). Since the binary components of GW190521 are unlikely both to have formed directly from stellar collapse, it is possible that they may be of PBH origin (Bird et al. 2016); however, theoretical expectations of the mass distribution and merger rate of PBH binaries have large uncertainties (e.g., Byrnes et al. 2018), so we do not attempt to quantify such scenarios. Some theories of PBH formation predict predominantly small component spins χ ≪ 1 (Chiba & Yokoyama 2017), which are disfavored by our parameter estimates for GW190521 (Section 2.2 and Figure 2).

Cosmic string cusps and kinks (e.g., Damour & Vilenkin 2000) may yield short-duration transient signals (bursts) with support at low frequency (e.g., Divakarla et al. 2019, Section IV D). There has been no previous detection of a cosmic string GW burst (Aasi et al. 2014b; Abbott et al. 2018b), and bounds on cosmic string model parameters derived from the overall contribution to the GW stochastic background are generally more stringent by orders of magnitude than bounds from direct burst searches (Abbott et al. 2018b, 2019a). Thus, detection of a burst signal in current data is a priori unlikely; however, we consider this possibility for completeness. Here we estimate the likelihood for the data to be produced by a cosmic string signal and compare this with binary BH merger models as described in previous sections.

GW190521 is identified by the cosmic string matched filter search pipeline (Aasi et al. 2014b; Abbott et al. 2018b); however, the maximum S/Ns in this search (∼6 and ∼8 in LIGO Hanford and Livingston, respectively) are much lower than for modeled BBH search templates, best-fit binary merger waveform models, or unmodeled reconstructions, suggesting that the data strongly prefer a binary merger model to a cosmic string or cusp. In addition, the signal appears inconsistent with the cosmic string template, as evidenced by large values of the χ² test statistic, over 7 (3) in LIGO Livingston (Hanford), while the expected value for a cosmic string signal is unity. Figure 13 illustrates the mismatch between the best-matching cusp waveform and the Livingston data.

Zoom In Zoom Out Reset image size

As a refinement of this analysis accounting for parameter uncertainties and phase-space volumes for different models, we compare the Bayes factors for a BBH waveform model to cosmic string cusp and kink models, respectively. Since the true rates of binary mergers in the mass range of GW190521 and of cosmic string transient signals are both unknown, we do not attempt to quantify prior odds for the two hypotheses. We use methods and models, including priors on cosmic string signal parameters, detailed in Divakarla et al. (2019). For the BBH model, we use the NRSur PHM waveform model with priors on component masses, spins, and distance matching those in Section 2. Using the dynesty nested sampling algorithm (Speagle 2020) within Bilby (Ashton et al. 2019; Romero-Shaw et al. 2020), we find Bayes factors ${\mathrm{log}}_{10}{B}_{\mathrm{cusp}}^{\mathrm{BBH}}\,=29.8$ and ${\mathrm{log}}_{10}{B}_{\mathrm{kink}}^{\mathrm{BBH}}=29.7$, strongly favoring the BBH hypothesis; thus, it is highly unlikely that this event was of cosmic string origin.

CCSN signals modeled in multidimensional simulations (e.g., Mezzacappa et al. 2020; O'Connor & Couch 2018; Andresen et al. 2019; Powell & Müller 2019) are of long duration (up to 1.5 s), may be broadband (from few Hz to around 1–2 kHz), and do not have morphology similar to GW190521. Such signals are expected to be detectable for sources within a distance of a few tens of kiloparsecs from Earth (Abbott et al. 2020d). We expect a CCSN signal to have electromagnetic and neutrino counterparts; however, no counterpart has been observed within a few hundred kiloparsecs (IceCube Collaboration 2019). GWs from extreme emission CCSNe, produced by highly deformed, rapidly rotating or fragmented cores, could be observed at a range of a few tens of megaparsecs (Abbott et al. 2020d). At that distance we do not expect to detect neutrinos; however, a detection of GWs produced by these extreme emission scenarios is overall very unlikely, and an electromagnetic counterpart would still be expected. It is thus implausible that the GW190521 signal was produced by CCSNe.