The gravitational-wave (GW) signal GW190521 Abbott et al. 2019a ) was observed by the Advanced LIGO (Aasi et al. 2015) and Advanced Virgo (Acernese et al. 2015) detectors during their third observing run (O3). The event was found with four different search pipelines, both at low-latency and offline with improved background estimation; an offline search sensitive to generic transients found GW190521 with a three-detector network signal-to-noise ratio of 14.7, and an estimated false-alarm rate of 1 in 4900 yr . Another candidate GW signal was reported later on the same day (Abbott et al. 2019b ). The source of GW190521 is consistent with being a high mass binary black hole (BBH) system. The final merger product of GW190521, with an estimated mass of 142 +28 −16 M (all values quoted as medians with symmetric 90% credible interval), is the first strong observational evidence for an intermediate-mass black hole (IMBH) in the mass range 10 2 − 10 3 M , under the assumption of a quasi-circular BBH coalescence. This merger of two high mass black holes (primary mass 85 +21 −14 M , secondary mass 66 +17 −18 M ) is also exceptional as the first observation of a black hole (BH) that lies with high confidence in the mass gap predicted by pair-instability (PI) supernova theory (Woosley 2017) ; the probability that the primary mass is below 65 M is 0.3%. This high component mass represents a challenge for current astrophysical formation scenarios.
The very short duration (approximately 0.1 s) and bandwidth (around 4 cycles in the frequency band 30-80 Hz) of GW190521 means that the interpretation of the source as being a quasi-circular compact binary coalescence consisting of inspiral, merger and ringdown phases is not certain. Under that interpretation we find that the observed signal, including its frequency evolution, is entirely consistent both with three different waveform models derived from analytical and/or numerical solutions of general relativity (GR), and with direct comparisons to numerical relativity (NR) solutions. Therefore, most of the discussion in this paper and in proceeds under that assumption, including our inferences on the inferred masses and spins, and on the effect of including higher-order multipoles in the waveform models. However, as discussed below, other interpretations are possible, adding to the exceptional nature of this event.
Searches for IMBH binaries with total mass >100 M and primary mass 500 M were carried out in data from Initial LIGO and Virgo (Abadie et al. 2012; Aasi et al. 2014a ) and from the first and second observing runs of the Advanced detector era, O1 and O2 (Abbott et al. 2017a (Abbott et al. , 2019c . However, no significant candidates were identified: see Udall et al. (2020) for further discussion. The most stringent upper limit on the local IMBH merger rate from O1 and O2 is 0.20 Gpc −3 yr −1 (in comoving units, 90% confidence level), for binaries with equal component masses m 1 = m 2 = 100 M (Abbott et al. 2019c ). Other groups have also searched LIGO-Virgo open data (Abbott et al. 2018a (Abbott et al. , 2019d for possible IMBH events (Zackay et al. 2019; Nitz et al. 2019) . The O3 run started in April 2019 with significantly increased sensitivities for all three Advanced detectors compared to O1 and O2 (Tse et al. 2019; Acernese et al. 2019 ); here we consider the implications of GW190521, detected in the first half of the run, O3a (1 April through 1 October, 2019).
Observational evidence for IMBHs, usually defined as BHs with mass 10 2 − 10 5 M (see e.g., van der Marel 2004; Miller & Colbert 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. 2020) . 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; Lodato & Natarajan 2006; Begelman et al. 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; Portegies Zwart et al. 2004; Gurkan 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; Lützgendorf et al. 2011; Miller-Jones et al. 2012; Nyland et al. 2012; Strader et al. 2012; Lützgendorf et al. 2013; Lanzoni et al. 2013; Perera et al. 2017; Kızıltan et al. 2017; Lin et al. 2018; Tremou et al. 2018; Baumgardt et al. 2019; Zocchi et al. 2019; Mann et al. 2019) .
Several ultra-luminous X-ray sources, defined as those with a total luminosity, assumed isotropic, of 10 39 erg s −1 , have been studied as IMBH candidates (Kaaret et al. 2001; Matsumoto et al. 2001; Strohmayer & Mushotzky 2003; van der Marel 2004; Miller & Colbert 2004; Feng & Soria 2011; Sutton et al. 2012; Mezcua et al. 2015; Wang et al. 2015; Mezcua 2017; 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 Godet et al. 2009; Servillat et al. 2011; Webb et al. 2012; Cseh et al. 2015; Soria et al. 2012) , pointing to an IMBH mass ∼0.3 − 30 × 10 4 M .
Several IMBH candidates lie at the centers of dwarf galaxies and are associated with low-luminosity activegalactic nuclei (AGNs; Filippenko & Sargent 1989; Filippenko & Ho 2003; Barth et al. 2004; Greene & Ho 2004; Barth et al. 2005; Greene & Ho 2007; Seth et al. 2010; Dong et al. 2012; Reines et al. 2013; Baldassare et al. 2015; den Brok et al. 2015; Baldassare et al. 2016; Mezcua et al. 2016; Baldassare et al. 2017; Mezcua et al. 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 ∼10 2 −10 3 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 pair instability (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 or with direct collapse, leaving a compact object less massive than expected in the absence of PPI (Woosley 2017 (Woosley , 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; Giacobbo et al. 2018; Woosley 2019) . The boundaries of the mass gap are highly uncertain because they depend on stellar evolution and on our understanding of core-collapse, PPI and PI supernovae (Stevenson et al. 2019; Marchant et al. 2019; Farmer 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 , 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 quasicircular 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, and 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 implica-tions 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 core-collapse supernova as possible sources due to mismatch with the GW data. Finally, we summarize our observations and consider future prospects in Section 7.
Under the assumption that GW190521 is a quasicircular BBH coalescence, the intrinsic parameters of the source are fully described by the masses, m 1 and m 2 , and the spin vectors, S 1 and S 2 , of the two BHs. We use the convention that m 1 ≥ m 2 and the massratio q = m 2 /m 1 ≤ 1. The dimensionless spin magnitudes, χ i = |c S i /(Gm 2 i )|, are assumed to be constant throughout the inspiral, while the spin orientations relative to the orbital angular momentum axis, θ LSi = arccos(L·Ŝ i ) ∈ [0, 180 • ], evolve over the duration of the signal. The spin orientations must therefore be parameterized at some fiducial time which, for this work, is when the signal has a GW frequency f 0 = 11 Hz. The remnant BH produced from post merger is described by its mass M f and spin magnitude χ f . We also estimate the BH recoil velocity v f = |v f | relative to the center of mass of the binary.
As discussed in Section 2.4 below, the source is estimated to be cosmologically distant. The masses measured in the frame of LIGO and Virgo detectors are therefore redshifted by a factor (1 + z) and are denoted with a superscript det so that m det = (1 + z) m, where m is the source frame mass. GWs directly encode the luminosity distance to the source D L , which in turn depends on the inclination angle of the binary orbit with respect to the line of sight (Section 2.4). To make inferences about the source frame masses we must therefore convert the distance measurement to a redshift. The statistical uncertainty associated with estimation of the source frame masses is increased relative to that of the detector frame masses due to these dependencies. From the inferred posterior distribution of D L we compute redshift assuming a Planck 2015 ΛCDM cosmology with Hubble parameter H 0 = 67.9 km s −1 Mpc −1 (Ade et al. 2016 ) (and we address the effect of taking a larger value of the Hubble parameter in Section 2.2). Unless stated otherwise, mass measurements are quoted in the source frame of the binary. 
To infer the source properties of GW190521, we analyzed 8 seconds of data in the LIGO and Virgo detectors around the time of the detection. The data are downsampled from 16384 Hz to 1024 Hz, as we expect no signal power above several hundred Hz due 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 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 . 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 Damour et al. 2001; Blanchet et al. 2005; Arun et al. 2009; Blanchet 2014 ) with tuning to NR solutions (Khan et al. 2020) ; SEOBNRv4PHM (SEOBNR PHM), an inspiralmerger-ringdown waveform model that is based on the effective-one-body formalism (Buonanno & Damour 1999 and calibrated to NR (Ossokine et al. 2020) .
These three waveform models employ different approaches to reproduce the predictions from analytical and numerical relativity; we expect to see differences in the parameter estimation from these three models due to 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 is 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 2 /m 1 ≥ 1/4. It includes all (l, |m|)-multipoles of the gravitational radiation (Blanchet et al. 1996; Kidder 2008; Blanchet et al. 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 . 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 doublespin model of spin-orbit precession during the inspiral (Chatziioannou et al. 2017) . After this precession "twisting", all l ≤ 4 modes will be non-zero in the inertial frame. The non-precessing, 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 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 sphericalspheroidal mode-mixing effects which can significantly impact some of the higher multipoles (Kelly & Baker 2013) .
SEOBNR PHM is based on the dynamics of spinning, non-precessing 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 non-precessing (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 non-precessing 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 Ossokine et al. 2020) . After the precession twisting, all l ≤ 5 modes will be non-zero 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. (2013) ; Lovelace et al. (2016) . The agreement between different NR codes is sufficiently good Pürrer & Haster (2020) to avoid systematic biases at the SNR 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 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 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 2 L (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 (right ascension and declination) 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 ini- tial 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) . gate 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 numerical relativity (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.
Masses -The estimated mass posterior distributions is shown in Fig. 1 (left) for the three GW signal models. The primary BH mass of GW190521 is m 1 = 85 +21 −14 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 2 = 66 +17 −18 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. 2019e) ; the secondary BH of GW190521 is also more massive than any primary BH in GWTC-1.
These source frame masses have been redshiftcorrected, as discussed above, using a value of the Hub-ble parameter H 0 = 67.9 from Planck 2015. However, recent measurements of H 0 using nearby Cepheid distance standards obtain a precise value of H 0 = 74.03 ± 1.42 km s −1 Mpc −1 (Riess et al. 2019) , 9% higher than the Planck value. Using this 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 1 and m 2 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 NR-Sur PHM and Phenom PHM models are 0.9% and 2.3%, respectively.
The probability that at least one of the black holes in GW190521 is in the range 65 − 120 M is 99.0%, using the NRSur PHM model. The corresponding probabilities using the Phenom PHM and SEOBNR PHM models are 98.0% and 90.2%, respectively.
We measure the total binary mass of GW190521 to be M = 150 +29 −17 M making it the highest-mass binary observed via GWs to date. The binary chirp mass is M = 64 +13 −8 M , a factor ∼2 times higher than the first BBH detection, GW150914 . GW190521 is consistent with a nearly equal mass binary with mass ratio q = m 2 /m 1 = 0.79 +0.19 −0.29 (90% credible interval).
In the detector frame, the measured masses are m det
−14 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 4 cycles (2 orbits) in the frequency band 30 − 80 Hz ), 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 log 10 -Bayes factor of +1.06 ± 0.06 for the NRSur PHM model allowing generic BH spins vs. limiting the effects of spin to the aligned components.
In the disk plots of Fig. 2 , we show constraints on the spins of the component BHs of GW190521 in terms of their dimensionless magnitudes χ 1 and χ 2 and polar angles (tilts) with respect to the orbital angular momentum, θ LS1 and θ LS2 , 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 − 0.9. The constraints on the tilt angles of these spins are also relatively broad.
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, χ eff = (m 1 χ 1 cos θ 1 + m 2 χ 2 cos θ 2 )/(m 1 + m 2 ) ∈ [−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 axisL 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 Fig. 1 (right) 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 towards 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. Fig. 3 ). The inferred mass of the remnant BH provides observational evidence for an IMBH of 100 M . The remnant BH mass is 7.6 +2.2 −1.9 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.7 −0.9 × 10 56 erg s −1 .
The remnant BH of GW190521 has a dimensionless spin parameter χ f = 0.72 +0.09 −0.12 . Within the 90% credible interval this is consistent with the BBH merger remnant spins reported in GWTC-1 (Abbott et al. 2019e ). 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 −1 for precessing systems 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 the posterior is weighted towards 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 ).
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 May 21, 2019 03:08:32 UTC (Abbott et al. 2019a ). The source was localized within a 90% credible area 1163 deg 2 (see blue curve in Fig. 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 2 (orange curve in Fig. 5 ). Here we report our latest constraints on the sky position of GW190521 found using LALInference with the NR-Sur PHM model (green curve in Fig. 5 ). The source is now localized within 774 deg 2 .
Distance and Inclination -The inferred luminosity distance D L reported at low latency (Abbott et al. 2019a ,f) 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 which 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 in-clination 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 towards 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
The masses of the component 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 close to face-on binary, 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.
The GWs emitted during a BBH coalescence can be decomposed as spherical harmonic multipoles and |m| ≤ .
The quadrupole ( , |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 for instance Calderón Bustillo et al. 2016; Pan et al. 2011 ). Higher-order multipoles have been found to be important for high mass BBH mergers in advanced detector observations (Varma et al. 2014; Graff et al. 2015; Calderón Bustillo et al. 2017; Calderón Bustillo et al. 2018; Chatziioannou et al. 2019; Pan et al. 2011; Mehta et al. 2017; 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 ) and GW190814 (Abbott et al. 2020b) . To quantify the impact and evidence of higher-order multipoles on GW190521, we computed the posterior distribution of our preferred model, NR-Sur PHM (Varma et al. 2019) , using the RIFT pipeline
. Sky maps (source location 90% credible areas) for GW190521, as seen from the north (left) and south (right) celestial poles. The blue (dashed) and the orange (solid thin) curves show two low-latency sky maps from the BAYESTAR pipeline and the LALInference pipeline using the SEOBNRv4 ROM waveform model; neither incorporates higher-order multipoles. The green (solid thick) curve, reported here for the first time, was obtained from full parameter estimation with the NRSur HM model. 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 .
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 Fig. 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 1 = 137 M at 90% credible interval. As we include higher-order multipoles up to = 4, the masses are better constrained and we find that m 1 cannot be greater than 103 M at 90% credible interval. We find that the higher-order multipoles have marginal impact on constraining the effective precession spin parameter χ p .
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 10 -Bayes factor for our preferred model with and without higher-order multipoles is −0.38, implying that the data marginally disfavors 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; Healy et al. 2019; Boyle 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 , the likelihood of the data was calculated for each waveform derived from NR simulations. Using the RIFT pipeline ) 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 Fig. 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 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.
GW190521 presents an opportunity for strong-field tests of general relativity in a previously unexplored region of parameter space (Abbott et al. 2016; Abbott et al. 2019g) . We first test the consistency of the data with the waveform models employed for parameter estimation by searching for unmodelled residual power; we then study in detail the properties of the ringdown phase using fewer assumptions on the progenitors of the remnant BH, and demonstrate consistency with predictions from the full waveform models.
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. 2016 . The residual data is analyzed by the BayesWave algorithm which simultaneously fits models for the noise PSD and short (< 1 s) transient signals coherent across the detector network . 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, SNR 90 , is the upper 90% credible bound on the posterior distribution function of the recovered signal-tonoise ratio of the wavelet reconstruction (Abbott et al. 2016) . The significance of the recovered SNR 90 value is empirically measured by repeating the analysis on 200 randomly chosen off-source times drawn from 4096 seconds of detector data surrounding the event time, which serves as the background estimate for SNR 90 . For GW190521, we find SNR 90 ∼ 5.97, which is consistent with expectations for typical LIGO-Virgo noise, resulting in a p-value of ∼0.4 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 general relativity, 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. Thse are extensions to similar tests applied to previous LIGO-Virgo BBH events (Abbott et al. 2016; Abbott et al. 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 Berti et al. ( , 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 which makes no assumptions on the process leading to the signal formation, and also waveform templates specifically modelling 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 modelling 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 0 (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 NR-Sur PHM model and yields a GPS time for LIGO Hanford t H peak = 1242442967.4306 +0.0067 −0.0106 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 Fig. 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. (2016) : we choose start times ∆t 0 ≡ t 0 − t H peak 6.4, 12.7, 19 .1 ms after the reference time. In units of the redshifted remnant mass M det f = (1 + z)M f = 258.3 M (median estimate from NRSur PHM) these times correspond to ∆t 0 = 5, 10, 15 GM det f /c 3 . Truncating the data at later times yields uninformative posteriors due to the decrease in SNR. Values of ∆t 0 ≥ 10 GM det f /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. The filled contours show the measurement using the fundamental Kerr = 2, m = 2, n = 0 multipole (black); the = 2, m = 2 Kerr model including overtones up to n = 2 (gray); and the fundamental higher mode model (Kerr HMs, light gray) described in the text. These are compared with redshifted remnant mass and spin obtained using the three waveform models stated in Section 2.1 and Fig. 3 . Contours enclose 90% of the posterior distribution, and the 1D histogram shows the 90% credible regions only for the ringdown models.
At a start time of ∆t 0 = 12.7 ms (corresponding to 10 GM det f /c 3 ), we measure a ringdown frequency f = 66.0 +4.0 −3.0 Hz and damping time τ = 19.0 +9.0 −7.0 ms. As the start time increases from 6.4 to 19.1 ms, we observe convergence towards the predicted value of the leastdamped mode frequency, together with a broadening of the posterior due to the decreasing SNR 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) ; Okounkova (2020) for discussions concerning the interpretation of a linearised 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 Fig. 9 , assuming a uniform prior over these two parameters. We 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 sphericalspheroidal 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 numerical relativity simulations of non-precessing BBH mergers. With this model we measure M det f = 296.0 +52.0 −58.0 M , χ f = 0.79 +0.12 −0.3 taking ∆t 0 = 12.7 ms; the full probability distribution is shown in the right panel of Fig. 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 analy-ses. 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 favour of the presence of higher multipoles or overtones.
A parametrized test of gravitational waveform generation (Li et al. 2012a,b; Agathos et al. 2014; Blanchet & Sathyaprakash 1995; Mishra et al. 2010 ) based on the Phenom PHM waveform model, also does not reveal inconsistencies with GR predictions. Full details will be provided in an upcoming paper.
We estimate the rate of mergers similar to this source, assuming a constant rate per comoving volume-time element. We proceed similarly to : in the absence of a parameterized population model for such sources, we assume a population of mergers whose intrinsic parameters (component masses and spins) are identical to the detected event up to measurement errors (Kim et al. 2003) . We estimate the sensitivity of the LIGO-Virgo detector network by adding simulated signals to data from the O1, O2 and O3a observing runs and recovering them with the Coherent WaveBurst (cWB) weakly-modeled transient detection pipeline (Klimenko et al. 2016) , optimized for sensitivity to IMBH mergers (Abadie et al. 2012) , which identified GW190521 with the highest significance . As in Abbott et al. (2017c) and Abbott et al. (2016b) , we consider a simulated signal to be detected if recovered with an estimated false alarm rate of 1 per 100 years or less. Simulated signals are prepared by drawing source parameter samples from a posterior distribution inferred using the NRSur PHM waveform (Section 2). The simulations' component masses and spins are taken directly from the posterior samples, whereas their lineof-sight direction and orbital axis direction are randomized, and their luminosity distances are distributed uniformly over comoving volume and time.
The time-and angle-averaged sensitive luminosity distances of the cWB-IMBHB search for mergers similar to GW190521 in O1, O2 and O3a data are 1.1 Gpc, 1.2 Gpc and 1.7 Gpc, respectively, showing a substantial gain in sensitivity over successive runs. The combined searched time-volume over O1, O2 and O3a data is 9.1 Gpc 3 yr. Thus, taking a Jeffreys prior on the rate R as p(R) ∝ R −0.5 , and with 1 event detected above threshold, we obtain an estimate R = 0.13 +0.30 −0.11 Gpc −3 yr −1 . This rate is below previous upper limits obtained by various methods: (mass-and spin-dependent) up-per limits inferred from LIGO-Virgo searches of O1-O2 data, 0.20 Gpc −3 yr −1 or greater (Abbott et al. 2019c ); limits obtained in Chandra et al. (2020) from numerical simulations of BBH signals with generic precessing spins added to O1 data, 0.36 Gpc −3 yr −1 or greater. Our rate is also well below model-dependent limits of O(1) Gpc −3 yr −1 obtained in Fishbach et al. (2020) for possible BBH populations with 45 < m 1 /M < 150. Our estimate is for a population sharing the properties of GW190521; as we obtain more observations of high mass BBH systems we expect to better constrain the distribution of their masses and spins, and thus refine the population rate estimate. The merger rate obtained here may be compared to expectations from possible formation channels, as in the following section.
Both the primary component mass and remnant mass of this system are higher than the previously most massive BBH detected by LSC and Virgo, GW170729 (Abbott et al. 2019e). Other candidate events, dubbed GW170817A and GW151205, have been identified in O1 and O2 data by other groups (Zackay et al. 2019; Nitz et al. 2019) , with component masses higher than GW170729 if of astrophysical origin. Here we do not attempt to assess possible astrophysical implications of such additional events.
Our analysis of the BBH population detected in the O1 and O2 runs using parameterized models of the mass distribution indicates that 99% of merging BH primaries have masses below 45 M (Abbott et al. 2019h; . As the primary mass of GW190521 is well above this value, the system is well within the highest mass 1% of the population inferred from O1 and O2 observations. We obtain confirmation of the unexpected nature of GW190521 by generating synthetic catalogs of 25 (50) BBH detections to represent expectations for BBH detections in the first half of O3. The masses of synthetic detections are obtained via draws from the posterior of the O1-O2 population described by the most general model, Model C, of Abbott et al. (2019h) , after applying selection effects. We then extract the highest-mass primary BH m 1,max and compare to the primary mass of GW190521. For the 25-event (50event) synthetic catalogs, 0.6% (1.3%) of m 1,max values lie above 85 and 3.0% (5.3%) lie above 71 (the posterior median and 5th percentile estimated using the NR-Sur PHM model). Even when accounting for statistical uncertainties in the mass estimation (Fishbach et al. 2020) , the primary mass of GW190521 is in tension with the population inferred from O1 and O2. Figure 10 shows the mass of GW190521 in comparison with the masses of all the O1 and O2 BBHs (left-hand panel) and with current theoretical knowledge about the PI mass gap (right-hand panel). A set of evolutionary astrophysical models (Spera & Mapelli 2017 ) relating the progenitor mass and compact object is also shown for reference. The astrophysical models are subject to several uncertainties in stellar evolution and corecollapse supernovae, and only serve as representative examples. This figure is an update of Fig. 5 of Abbott et al. (2019h) to include the masses of GW150521.1 and the most recent uncertainty estimates on the PI mass gap, which we review below in Section 5.1.
Several mechanisms could fill the PI gap. In the following sections, we discuss second-generation (2g) BHs, stellar mergers in young star clusters and BH mergers in AGN disks. Second-generation BHs, i.e. BHs born from the merger of two BHs, can have mass in the PI gap (Miller & Hamilton 2002; Gerosa & Berti 2017; Rodriguez et al. 2019 ). If they are in a dense stellar environment and are not ejected by the gravitational recoil, they have a chance to form a new binary system with another BH Gerosa & Berti 2019 ). Alternatively, a merger between an evolved star (with a well developed helium or carbon-oxygen core) and a main-sequence companion might trigger the formation of a giant star with an over-sized envelope with respect to the core. If this star collapses into a BH before its helium core enters the PI range, then it can give birth to a BH in the PI gap (Spera et al. 2019; Di Carlo et al. 2019) . If this BH is inside a dense stellar environment, it has a chance to capture a companion by dynamical exchange. BHs in AGN disks might pair with other BHs and are expected to merge efficiently due to gas torques, producing 2g BHs (e.g., McKernan et al. 2018) . These BHs in AGN disks might even grow by gas accretion (e.g., McKernan et al. 2012) .
Finally, van Son et al. (2020) investigated the possibility that BHs in isolated binaries might have mass in the PI gap by super-Eddington accretion of the companion onto the primary BH. Even under the most extreme assumptions, they find that at most ∼ 2% of all merging BBHs in isolated binaries contains a BH with a mass in the pair-instability mass gap, and find no merging BBH with a total mass exceeding 100 M .
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 Woosley , 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 which combine with PI and PPI to shape the mass spectrum of BHs (e.g., Belczynski et al. 2016; Stevenson et al. 2019; Belczynski et al. 2020) . 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 modelling 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 12 C(α, γ) 16 O reaction ; see also Takahashi et al. 2018 ). This result is particularly important if we consider that Farmer et al. 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 −4 ) 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 secondgeneration 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. 2019e ) 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 sin-gle 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 Fig. 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 −1 (Holley-Bockelmann et al. 2008; O'Leary et al. 2016; Antonini & Rasio 2016; Gerosa & Berti 2019) . In this scenario, the companion might be a 1g BH or, less likely, a 2g BH .
The relativistic kick at birth strongly depends on BH spin: maximally spinning BHs with spins in the orbital plane and counter-aligned receive the largest kicks (e.g., Campanelli et al. 2007a,b) . 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 . For these reasons, hierarchical mergers are expected to be more likely in nuclear star clusters (e.g., Antonini & Rasio 2016; 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 (Gerosa & Berti 2019; Rodriguez et al. 2019; Chatziioannou 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 . 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 BH where one component was above the PI mass gap, though this is disfavored by the preference for near-equal masses (see Fig. 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; ). 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 (Gerosa & Berti 2019; Rodriguez et al. 2019; Kimball et al. 2020a; Farr et al. 2018 ). 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 −1 , 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.
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 log 10 P (m1,χ1|xg+2g) P (m1,χ1|1g+1g) Figure 11 . Marginalized Bayes factor for a BBH to be a second-generation merger (1g+2g or 2g+2g) as opposed to a first-generation BH merger, as a function of primary mass and spin, in the globular cluster analysis of Section 5.2.1 using a physically-motivated prior cutoff on 1g BH masses. The Bayes factor contours correspond to component masses and spins inferred using the NRSur PHM model for GW190521, and differ only slightly from those found using the Phenom PHM model. We show the 90% and 68% posterior credible regions for GW190521 as solid and dashed contours, respectively, for both the NRSur PHM and Phenom PHM models.
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 non-singular 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 ) and fallback dynamics 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 which 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 5 M with a Plummer radius of 1 pc (Harris 1996; Kremer et al. 2020) ; to quantify dependence on the escape velocity, we also consider a higher mass cluster of 10 8 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 BBH 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 sub-population of 1g BHs with zero spin, making up a fraction λ 0 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 vs 1g+1g origin as a function of primary component mass and spin, together with the 90% and With zero-spin channel Without zero-spin channel Figure 12 . Relative rates of 1g+2g and 2g+2g as compared to first-generation mergers, in globular cluster models with (blue) and without (orange) a zero-spin stellar BH population (see Section 5.2.1), using GW190521 source parameters derived from NRSur PHM. In the model with zero-spin population, we also plot the fraction λ0 of 1g+1g binary components belonging to this population.
68% credible regions for GW190521. Overall Bayes factors including dependence on both component masses and spins generally favor 2g+2g origin; the parameterestimation likelihood strongly favors near-equal component masses, thus 1g+2g origin implying significantly unequal masses is disfavored by the data. Figure 12 (orange contours) shows the relative rates of 1g+2g and 2g+2g mergers for the model without a 1g zero-spin sub-population.
We find log 10 R 1g+2g /R 1g+1g and log 10 R 2g+2g /R 1g+1g to be −2.6 +0.9 −1.1 and −5.5 +1.8 −2.2 , 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 BBH with large component spins are often subject to kicks which 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.
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 Fig. 12 (blue contours): log 10 R 1g+2g /R 1g+1g and log 10 R 2g+2g /R 1g+1g become −2.0 +1.0 −1.2 and −4.4 +2.0 −2.4 , 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, both for the NRSur PHM and Phenom PHM parameter estimates. In this case, the estimated fraction λ 0 of merging 1g BHs from the zero-spin channel is in the range ∼0.01-0.47.
We also consider a higher cluster mass of 10 8 M which may be representative of nuclear cluster enviroments: 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 modelling of BH mergers in AGN 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 over-sized 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 (Spera et al. 2019; Di Carlo et al. 2019; Bouffanais 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 over-sized 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 −5 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 due 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 ef-ficiently torqued by gas drag till they align with the AGN disk (Bartos et al. 2017; McKernan et al. 2018) . Once in the disk, BHs can accrete gas, and are expected to acquire companions and efficiently merge with them as an effect of gas torques (McKernan et al. 2012 (McKernan et al. , 2014 Bellovary et al. 2016; Bartos et al. 2017; McKernan et al. 2020) . 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 km s −1 , given the proximity to a SMBH). Hence, AGN disks might easily harbor 2g BHs (Yang et al. 2019b) .
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. 2019a) , or even reminiscent of a broken power law (McKernan et al. 2018; Yang et al. 2019b ). 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. 2019b) .
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 ultra-violet flare; McKernan et al. 2019) . No prompt electromagnetic counterpart has been reported for GW190521 (e.g., Lipunov et al. 2019; Doyle et al. 2019a; Casentini et al. 2019; Doyle et al. 2019b; Negoro et al. 2019; Watson et al. 2019; Turpin et al. 2019a,b; Barthelmy 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. (2019b) 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 ∼1 − 10 Gpc −3 yr −1 (e.g., Bartos et al. 2017; Stone et al. 2017; Yang et al. 2019a; Sedda 2020 ), while McKernan et al. (2018 and Tagawa et al. (2020) try to quantify additional uncertainties and end up with a merger rate ∼10 −3 − 10 4 Gpc −3 yr −1 and ∼0.02 − 60 Gpc −3 yr −1 , respectively.
The short duration and bandwidth of the signal observed in the data open the possibility that the source may not be uniquely explained as an unusually massive, quasi-circular BBH merger formed via astrophysical processes. Furthermore, given a transient GW signal for which a full inspiral-merger-ringdown morphology is not directly evident in the data, a natural question is whether the signal may be consistent with other types of modelled GW source. We consider several possible alternative scenarios below; all are disfavored either by the data, or by low prior probability of the alternative hypothesis, or by both.
The waveform observed in the data may also be consistent with the merger of a BBH with non-zero 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 ). 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 Samsing & Ramirez-Ruiz 2017; Samsing 2018; Rodriguez et al. 2018; Samsing & D'Orazio 2018; Zevin et al. 2019; Fragione et al. 2020) . Head-on collisions are even more rare due 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; Silsbee & Tremaine 2017; Antonini et al. 2017) or in dense star clusters (Wen 2003; Antonini & Perets 2012; 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 −3 yr −1 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 non-zero 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 (Ng et al. 2018; Broadhurst et al. 2018; Oguri 2018; Li et al. 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 one in every hundred to thousand detected events strongly lensed by individual galaxies (Holz & Wald 1998; Ng et al. 2018; Oguri 2018; Li et al. 2018 ) and a similarly low rate for lensing by galaxy clusters Oguri 2019) .
Under the strong lensing hypothesis H L , we assume 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 )
where ϑ are the apparent parameters of the waveform received at the detector, which differ from the sourceframe parameters due to effects of redshift and lensing. For the prior over component masses and redshift of the source p(m 1 , m 2 , z|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 −6 (D c (z)/Gpc) 3 (1.8 +1.7 −0.6 ). At these redshifts, the lensing optical depth is low: 9 +10 −3 × 10 −4 (5 +9 −3 × 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. 2019b) , was reported 4.6 hours 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 .
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 sub-threshold 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 Fig. 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. 2020, Section IV D) . There has been no previous detection of a cosmic string GW burst (Abbott et al. 2018b; Aasi et al. 2014b ) 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 . 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 (Abbott et al. 2018b; Aasi et al. 2014b ); however, the maximum signal-to-noise ratios in this search (∼6 and ∼8 in LIGO Hanford and Livingston, respectively) are much lower than for either modelled BBH search templates, or best-fit binary merger waveform models, or for unmodelled reconstructions, suggesting the data strongly prefers 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 χ 2 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.
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. (2020) . 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 log 10 B BBH cusp = 29.8 and log 10 B BBH kink = 29.7, strongly favoring the BBH hypothesis, thus it is highly unlikely that this event was of cosmic string origin. Bandpassed and whitened data Best matching CS template Figure 13 . The LIGO Livingston strain data time-series is shown in dark grey. The best-matching cosmic string cusp template is shown in red. Both are whitened using the detector's noise spectrum; the strain data are additionally bandpassed between 20 Hz and 128 Hz to suppress noise well above the frequency support of the template.
Core-collapse supernova (CCSN) signals modeled in multidimensional simulations (e.g., Andresen et al. 2019; O'Connor & Couch 2018; Powell & Müller 2019; Mezzacappa et al. 2020 ) are of long duration (up to 1.5 s), may be broad-band (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 few tens of kpc distance from Earth ). We expect a CCSN signal to have electromagnetic and neutrino counterparts; however no counterpart has been observed within a few hundred kpc (S190521g GCNs 2019). GWs from extreme emission CCSNe, produced by highly deformed, rapidly rotating or fragmented cores, could be observed at a few tens of Mpc range . 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 CCSN.
The GW signal GW190521 observed by the Advanced LIGO and Advanced Virgo detectors is consistent with emission from the coalescence of a high mass BBH system. Under that assumption, at least one of the component BHs in the binary has a mass in the PI mass gap with high probability, and the final merged BH has a mass in the IMBH range. There is no conclusive previous evidence from electromagnetic observa-tions for the existence of IMBHs in this mass range (∼100 − 1000 M ).
The signal is well described by waveforms derived from general relativity incorporating spin-induced orbital precession and higher-order multipoles, but neglecting orbital eccentricity. There is only weak evidence for the effects of orbital precession in the signal, and there is moderate but not conclusive evidence for high spin components in the orbital plane. The Bayes factor for the presence of higher-order multipoles is slightly negative, disfavoring their presence; however their inclusion enables more precise estimates of the source's component masses, distance and inclination to the line of sight.
Several tests demonstrate the consistency of the signal with GR predictions, including agreement between full inspiral-merger-ringdown waveforms derived from numerical relativity, and ringdown-only waveforms formed by summing the quasi-normal modes of the final Kerr BH. Because of the short duration and low bandwidth of the observed signal, possible sources other than a quasicircular, high mass BBH merger may be considered. Initial studies using waveforms that include orbital eccentricity suggest a partial degeneracy between eccentricity significantly greater than 0.1 and precession; however, significant eccentricity is disfavored by a low prior probability. The possibility that the signal is from a cosmic string cusp or kink is considered, but the signal is a poor match to waveform models from such sources. Finally, the possibility of the source being a core collapse supernova is considered; this is disfavored by the lack of any counterpart electromagnetic or neutrino signal for CCSN.
We discussed various scenarios for the formation and evolution of such a massive system. Uncertainties on the PI mass gap might justify the formation of BHs with mass > 65 M from stellar collapse. Alternatively, the primary BH might be the result of the merger of two smaller BHs (hierarchical scenario), or of two massive stars. The mild evidence for precession suggests a dynamical formation scenario, which predicts nearly isotropic spin orientation, such as the hierarchical BH merger or the outcome of a stellar merger in star clusters. The formation of GW190521 via isolated binary evolution appears disfavored.
Finally, we considered the possibilities that GW190521 is a strongly lensed signal from a lowermass, high redshift BBH merger, or that the system is a binary of primordial BHs (formed in the early Universe). The mass of the system, the mass ratio and the merger rate derived for this event are consistent with all of these scenarios, although the prior probability of strong lensing is low.
It is not possible to conclude at this time whether GW190521 represents the first of a new population of BBHs, or is merely at the high-mass end of the population of BBH systems already observed by LIGO and Virgo. A future publication will address this question in the context of the larger sample of BBH events observed in the first half of the LIGO-Virgo O3 run, O3a (2019 April 1 through October 1). The answer, whether positive or negative, has the potential to provide new insights into the population and evolution of the most massive stars. We look forward to observing more binary mergers with high mass to further inform our understanding of these phenomena.
In upcoming future observing runs (Abbott et al. 2020b) we expect the global advanced detector network's sensitivity to BBH mergers to increase significantly, with potentially several hundreds of detections per year, reaching out to redshifts of a few or more. We may thereby observe a large sample of events similar to GW190521, which will constitute a unique source of information on the binary formation environments and channels, and on the dynamics and evolution of massive stellar BHs over cosmic history. This will enable us to address questions such as the natal spin and mass distribution of BHs as a product of stellar collapse dynamics including PI; stellar cluster masses and dynamics; merger kicks due to GW emission; and tests of general relativity for highly spinning BHs.
The evidence for high mass BBH and IMBH binaries advances the science case for enhanced sensitivity at lower GW frequencies ( 5 Hz) . The proposed nextgeneration ground-based GW detectors, such as Einstein Telescope (Punturo et al. 2010) and Cosmic Explorer (Abbott et al. 2017d; Reitze et al. 2019) , can observe events similar to GW190521 up to redshift ∼20 (Gair et al. 2011; . Such potential observations at high redshift can probe the dynamics and evolution of massive stellar BHs over cosmic history. GW190521 also motivates the possibility for multiband GW observations Miller 2009; Sesana 2016; ) between ground-based detectors and upcoming space-based detectors like LISA (Amaro-Seoane et al. 2017). Such joint observations will provide a unique opportunity to test general relativity with highly spinning BHs (Vitale 2016) and probe binary formation environments and channels through measurements of spin evolution and eccentricity (Cutler et al. 2019 ).
The authors gratefully acknowledge the support of the United States National Science Foundation (NSF) for the construction and operation of the LIGO Laboratory and Advanced LIGO as well as the Science and Technology Facilities Council (STFC) of the United Kingdom, the Max-Planck-Society (MPS), and the State of Niedersachsen/Germany for support of the construction of Advanced LIGO and construction and operation of the GEO600 detector. Additional support for Advanced LIGO was provided by the Australian Research Council. The authors gratefully acknowledge the Italian Istituto Nazionale di Fisica Nucleare (INFN), the French Centre National de la Recherche Scientifique (CNRS) and the Netherlands Organization for Scientific Research, for the construction and operation of the Virgo detector and the creation and support of the EGO consortium. The authors also gratefully acknowledge research support from these agencies as well as by the We would like to thank all of the essential workers who put their health at risk during the COVID-19 pandemic, without whom we would not have been able to complete this work.
Strain data from the LIGO and Virgo detectors associated with GW190521 can be found in . This document has been assigned the LIGO document number LIGO-P2000021.