The first two observing runs (O1 and O2) with Advanced LIGO (Aasi et al. 2015) and Advanced Virgo (Acernese et al. 2015) opened up the field of gravitational-wave astrophysics with the detection of the first binary black hole (BBH) coalescence signal, GW150914 (Abbott et al. 2016a) . Another nine such events (Abbott et al. 2016b (Abbott et al. , 2019a were discovered by the LIGO Scientific and Virgo Collaborations (LVC) during this period, and additional events were reported by independent groups (Venumadhav et al. 2020; Zackay et al. 2019a,b; Nitz et al. 2020 ). The first binary neutron star (BNS) coalescence signal, GW170817, was discovered during the second of these observing campaigns (Abbott et al. 2017a (Abbott et al. , 2019b . It proved to be a multi-messenger source with emission across the electromagnetic spectrum (Abbott et al. 2017b) , with implications for the origin of short gamma-ray bursts (Abbott et al. 2017c) , the formation of heavy elements (Abbott et al. 2017d; Chornock et al. 2017; Tanvir et al. 2017; Rosswog et al. 2018; Kasliwal et al. 2019; Watson et al. 2019) , cosmology (Abbott et al. 2017e, 2019c and fundamental physics (Abbott et al. 2017c (Abbott et al. , 2019d .
The first six months of the third observing run (O3) were completed between April 1 and September 30, 2019. The LVC recently reported on the discovery of GW190425, the coalescence signal of what is most likely a BNS with unusually large chirp mass and total mass compared to the Galactic BNSs known from radio pulsar observations (Abbott et al. 2020a) . Another discovery from O3 is that of GW190412, the first BBH coalescence with an unequivocally unequal mass ratio q = m 2 /m 1 of 0.28 +0.12 −0.06 (all measurements are reported as symmetric 90% credible intervals around the median of the marginalized posterior distribution, unless otherwise specified). It is also the first event for which highermultipole gravitational radiation was detected with high significance (Abbott et al. 2020d ).
Here we report on another O3 detection, GW190814, the signal of a compact binary coalescence with the most unequal mass ratio yet measured with gravitational waves: q = 0.112 +0.008 −0.009 . The signal was first identified in data from two detectors, LIGO Livingston and Virgo, on 2019 August 14, 21:11:00 UTC. Subsequent analysis of data from the full three-detector network revealed a merger signal with signal-to-noise ratio (SNR) of 25.
The primary component of GW190814 is conclusively a black hole (BH) with mass m 1 = 23.2 +1.1 −1.0 M . Its dimensionless spin magnitude is constrained to χ 1 ≤ 0.07. The nature of the 2.59 +0.08 −0.09 M secondary component is unclear. The lack of measurable tidal deformations and the absence of an electromagnetic counterpart are consistent with either a neutron star (NS) or a BH given the event's asymmetric masses and distance of 241 +41 −45 Mpc. However, we show here that comparisons with the maximum NS mass predicted by studies of GW170817's remnant, by current knowledge of the NS equation of state, and by electromagnetic observations of NSs in binary systems indicate that the secondary is likely too heavy to be a NS. Either way, this is an unprecedented source because the secondary's well-constrained mass of 2.50-2.67 M makes it either the lightest BH or the heaviest NS ever observed in a double compact-object system.
As in the case of GW190412, we are able to measure the presence of higher multipoles in the gravitational radiation, and a set of tests of general relativity with the signal reveal no deviations from the theory. Treating this event as a new class of compact binary coalescences, we estimate a merger rate density of 1-23 Gpc −3 yr −1 for GW190814-like events. Forming coalescing compact binaries with this unusual combination of masses at such a rate challenges our current understanding of astrophysical models.
We report on the status of the detector network and the specifics of the detection in Sections 2 and 3. In Section 4, we estimate physical source properties with a set of waveform models, and we assess statistical and systematic uncertainties. Tests of general relativity are described in Section 5. In Section 6, we calculate the merger rate density and discuss implications for the nature of the secondary component, compact binary formation and cosmology. Section 7 summarizes our findings.
At the time of GW190814, LIGO Hanford, LIGO Livingston and Virgo were operating with typical O3 sensitivities (Abbott et al. 2020a) . Although LIGO Hanford was in a stable operating configuration at the time of GW190814, the detector was not in observing mode due to a routine procedure to minimize angular noise coupling to the strain measurement (Kasprzack & Yu 2017) . This same procedure took place at LIGO Hanford around the time of GW170608; we refer the reader to Abbott et al. (2017f) for details of this procedure. Within a 5 min window around GW190814, this procedure was not taking place, therefore LIGO Hanford data for GW190814 are usable in the nominal range (Chatterji et al. 2004 ) of data containing GW190814, observed by LIGO Hanford (top), LIGO Livingston (middle), and Virgo (bottom). Times are shown relative to 2019 August 14, 21:10:39 UTC. Each detector's data are whitened by their respective noise amplitude spectral density and a Q-transform is calculated. The colorbar displays the normalized energy reported by the Q-transform at each frequency. These plots are not used in our detection procedure and are for visualization purposes only.
of analyzed frequencies. A time-frequency representation (Chatterji et al. 2004 ) of the data from all three detectors around the time of the signal is shown in Figure 1. We used validation procedures similar to those used to vet previous gravitational-wave events (Abbott et al. 2016c (Abbott et al. , 2019a . Overall we found no evidence that instrumental or environmental disturbances (Effler et al. 2015) could account for GW190814. However, we did identify low-frequency transient noise due to scattered light at LIGO Livingston, a common source of noise in all three interferometers (Nuttall 2018) . Scattered light features in the strain data are produced when a small fraction of the main laser beam reflects off a moving surface and is phase modulated before recombining with the main beam. This recombination can result in excess noise with the morphology of arches in the time-frequency plane; the frequency of this noise is determined by the velocity of the moving surface (Accadia et al. 2010) . Thunderstorms near LIGO Livingston around the time of GW190814 resulted in acoustic noise coupling to the detector and caused features in the strain data associated with scattered light (Abbott et al. 2019a) . In this instance, this form of noise affects frequencies up to 30 Hz from roughly 22 s to 8 s before and 0.2 s to 1.5 s after the detected time of GW190814, as seen in the middle panel of Figure 1 . Since this noise could bias the estimation of GW190814's source parameters, we used a starting frequency of 30 Hz to analyse LIGO Livingston data. Virgo was operating nominally and there are no quality issues in the Virgo data.
The LIGO and Virgo detectors are calibrated by photon pressure from modulated auxiliary lasers inducing test-mass motion (Karki et al. 2016; Acernese et al. 2018; Viets et al. 2018) . Over the frequency range 20-2048 Hz, the maximum 1σ calibration uncertainties for strain data used in the analysis of GW190814 were 6% in amplitude and 4 deg in phase for LIGO data, and 5% in amplitude and 7 deg in phase for Virgo data. These calibration uncertainties are propagated into the parameter estimation reported in Section 4 via marginalization.
3. DETECTION 3.1. Low-latency Identification of a Candidate Event GW190814 was first identified on 2019 August 14, 21:11:00 UTC as a loud two-detector event in LIGO Livingston and Virgo data (SNR 21.4 and 4.3) by the low-latency GstLAL matched-filtering search pipeline for coalescing binaries (Cannon et al. 2012; Privitera et al. 2014; Messick et al. 2017; Sachdev et al. 2019; Hanna et al. 2020) . Matched-filtering searches use banks (Sathyaprakash & Dhurandhar 1991; Blanchet et al. 1995; Owen 1996; Owen & Sathyaprakash 1999; Damour et al. 2001; Blanchet et al. 2005; Cokelaer 2007; Harry et al. 2009; Brown et al. 2013; Ajith et al. 2014; Harry et al. 2014; Capano et al. 2016b; Roy et al. 2017 Indik et al. 2018 ) of modeled gravitational waveforms (Buonanno & Damour 1999; Arun et al. 2009; Blanchet 2014; Bohé et al. 2017; Pürrer 2016) as filter templates. A Notice was issued through NASA's Gamma-ray Coordinates Network (GCN) 20 min later (LIGO Scientific Collaboration, Virgo Collaboration 2019a) with a two-detector source local-ization computed using the rapid Bayesian algorithm BAYESTAR (Singer & Price 2016 ) that is shown in Figure 2 . The event was initially classified as "Mass-Gap" (Kapadia et al. 2020; LIGO Scientific Collaboration, Virgo Collaboration 2019b) , implying that at least one of the binary merger components was found to have a mass between 3-5 M in the low-latency analyses.
Other low-latency searches, including the matchedfiltering based MBTA (Adams et al. 2016) and Py-CBC (Usman et al. 2016; Nitz et al. 2017; Nitz et al. 2018 Nitz et al. , 2019 pipelines, could not detect the event at the time as its SNR in Virgo data was below their single-detector detection thresholds. Test versions of MBTA and the additional matched-filtering pipeline SPIIR (Hooper et al. 2012; Liu et al. 2012; Guo et al. 2018) operating with a lower SNR threshold also identified the event with consistent attributes.
Shortly thereafter, reanalyses including LIGO Hanford data were performed using GstLAL and PyCBC. A coincident gravitational-wave signal was identified in all three detectors by both searches, with SNR 21.6 in LIGO Livingston, 10.6 in LIGO Hanford, and 4.5 in Virgo data (as measured by GstLAL, consistent with SNRs reported by PyCBC). Results of these 3-detector analyses were reported in a GCN Circular within 2.3 hours of the time of the event (LIGO Scientific Collaboration, Virgo Collaboration 2019c,d), providing a 3-detector localization (Singer & Price 2016) constraining the distance to 220-330 Mpc and the sky area to 38 deg 2 at the 90% credible level. Another GCN Circular (LIGO Scientific Collaboration, Virgo Collaboration 2019e) sent 13.5 hours after the event updated the source localization to a distance of 215-320 Mpc, the sky area to 23 deg 2 , and the source classification to "NSBH" (Kapadia et al. 2020; LIGO Scientific Collaboration, Virgo Collaboration 2019b) , indicating that the secondary had a mass below 3 M . These updated sky localizations are also shown in Figure 2 . The two disjoint sky localizations arise because the low SNR in the Virgo detector (4.5) means that the data are consistent with two different signal arrival times in that detector.
Several external groups performed multi-messenger follow-up of the source with observations across the electromagnetic spectrum (e.g., Lipunov et al. 2019; Gomez et al. 2019; Antier et al. 2020; Andreoni et al. 2020; Dobie et al. 2019; Watson et al. 2020; Ackley et al. 2020; Vieira et al. 2020 ) and with neutrino observations (e.g., Ageron et al. 2019 ; The IceCube Collaboration 2019). No counterpart candidates were reported. The nondetection is consistent with the source's highly un-equal mass ratio and low primary spin (LIGO Scientific Collaboration, Virgo Collaboration 2019d,e; Fernández et al. 2020; Morgan et al. 2020) . Tentative constraints placed by multi-messenger studies on the properties of the system, such as the ejecta mass and maximum primary spin (Andreoni et al. 2020; Ackley et al. 2020; Kawaguchi et al. 2020; Coughlin et al. 2020) or the circum-merger density (Dobie et al. 2019 ) assuming a neutron-star-black-hole (NSBH) source, may need to be revisited in light of the updated source parameters we present in Sec. 4.1.
The significance of GW190814 was estimated by follow-up searches using improved calibration and refined data-quality information that are not available in low latency. They also used longer stretches of data for better precision (Abbott et al. 2016b,c) . With LIGO Hanford data being usable but not in nominal observing mode at the time of GW190814, we used only data from the LIGO Livingston and Virgo detectors for significance estimation. GW190814 was identified as a confident detection in analyses of detector data collected over the period from August 7 to August 15, 2019 by the two independent matched-filtering searches GstLAL and Py-CBC, with SNR values consistent with the low-latency analyses. The production version of PyCBC for O3 estimates significance only for events that are coincident in the LIGO Hanford and LIGO Livingston detectors, and therefore an extended version was used for GW190814 in order to enable the use of Virgo data in significance estimation.
GstLAL and PyCBC use different techniques for estimating the noise background and methods of ranking gravitational-wave candidates. Both use results from searches over non time-coincident data to improve their noise background estimation Messick et al. 2017; Usman et al. 2016) . Using data from the first six months of O3 and including all events during this period in the estimation of noise background, GstLAL estimated a false-alarm rate (FAR) of 1 in 1.3 × 10 3 yr for GW190814. Using data from the 8-day period surrounding GW190814 and including this and all quieter events during this period in noise background estimation, the extended PyCBC pipeline ) estimated a FAR for the event of 1 in 8.1 yr. The higher FAR estimate from PyCBC can be attributed to the event being identified by the pipeline as being quieter than multiple noise events in Virgo data. As PyCBC estimates background statistics using non-coincident data from both detectors, these louder noise events in Virgo data can form chance coincidences (Veitch et al. 2015) , is shown in green. The purple contour indicates the final sky localization as presented in this paper, which constrains the source to within 18.5 deg 2 at 90% probability.
with the signal in LIGO Livingston data and elevate the noise background estimate for coincident events, especially when considering shorter data periods. All estimated background events that were ranked higher than GW190814 by PyCBC were indeed confirmed to be coincidences of the candidate event itself in LIGO Livingston with random noise events in Virgo. The stated background estimates are therefore conservative (Capano et al. 2016a) . We also estimate the background excluding the candidate from the calculation, a procedure that yields a mean-unbiased estimation of the distribution of noise events (Capano et al. 2016a; Abbott et al. 2016d) . In this case, with GstLAL we found a FAR of < 1 in 10 5 yr while with PyCBC we found a FAR of < 1 in 4.2 × 10 4 yr. With both pipelines identifying GW190814 as more significant than any event in the background, the FARs assigned are upper bounds. When data from LIGO Hanford were included, GW190814 was also identified by the unmodelled coherent Wave Burst (cWB) search that targets generic gravitational-wave transients with increasing frequency over time without relying on waveform models (Klimenko et al. 2008 (Klimenko et al. , 2016 Abbott et al. 2016e ). We found a FAR of < 1 in 10 3 yr of observing time against the noise background from LIGO Hanford and LIGO Livingston data, consistent with the other searches. Table 1 . Source properties of GW190814: We report the median values along with the symmetric 90% credible intervals for the SEOBNRv4PHM (EOBNR PHM) and IMRPhenomPv3HM (Phenom PHM) waveform models. The primary spin magnitude and the effective precession is given as the 90% upper limit. The inclination angle is folded to [0, π/2]. The last column is the result of combining the posteriors of each model with equal weight. The sky location of GW190814 is shown in Figure 2 .
We infer the physical properties of GW190814 using a coherent Bayesian analysis of the data from LIGO Livingston, LIGO Hanford and Virgo following the methodology described in Appendix B of Abbott et al. (2019a) . Results presented here are obtained using 16 s of data around the time of detection. We use a low-frequency cutoff of 20 Hz for LIGO Hanford and Virgo and 30 Hz for LIGO Livingston for the likelihood evaluations, and we choose uninformative and wide priors, as defined in Appendix B.1 of Abbott et al. (2019a) . The LALInference stochastic sampling software (Veitch et al. 2015) is the primary tool used to sample the posterior distribution. A parallelized version of the parameter estimation software bilby (pbilby; Smith & Ashton 2019; Ashton et al. 2019) is used for computationally expensive signal models. The power spectral density used in the likelihood calculations is a fair draw estimate calculated with BayesWave .
This signal is analyzed under two different assumptions: that it represents a BBH, or that it represents a NSBH. For the BBH analyses, two different waveform families are used, one based on the effective-one-body approach (EOBNR; Bohé et al. 2017; Babak et al. 2017; Cotesta et al. 2018; Ossokine et al. 2020 ) and the other on a phenomenological approach (Phenom; Khan et al. 2016; Husa et al. 2016; London et al. 2018; Khan et al. 2019; Khan et al. 2020) .
For the NSBH analyses, we use BBH waveform models augmented with tidal effects (Matas et al. 2020; ). Systematic uncertainties due to waveform modeling are expected to be subdominant compared to statistical errors (Huang et al. 2020) . When sampling the parameter space with the SEOB-NRv4 ROM NRTidalv2 NSBH (Matas et al. 2020) and IMRPhenomNSBH waveform models, we obtained posterior distributions for the secondary component's tidal deformability Λ 2 that are uninformative relative to a uniform prior in Λ 2 ∈ [0, 3000]. The absence of a measurable tidal signature is consistent with the highly unequal mass ratio (Foucart et al. 2013; Kumar et al. 2017 ) and with the relatively large secondary mass (Flanagan & Hinderer 2008) . The large asymmetry in the masses implies that the binary will merge before the neutron star is tidally disrupted for any expected NS equation of state (Foucart et al. 2013) . Given that the signal carries no discernible information about matter effects, here we present quantitative results only from BBH waveform models.
Our primary analyses include the effect of subdominant multipole moments in precessing waveform template models (PHM): IMRPhenomPv3HM (Phenom PHM; Khan et al. 2019; Khan et al. 2020 ) from the phenomenological family and SEOBNRv4PHM (EOBNR PHM; Babak et al. 2017; Ossokine et al. 2020 ) from the EOBNR family. 1 Analyses that assume the spins are aligned with the orbital angular momentum were also performed, either including (Phenom/EOBNR HM) or excluding (Phenom/EOBNR) the effect of subdominant multipole moments.
From the ∼ 300 observed cycles above 20 Hz, we are able to tightly constrain the source properties of GW190814. Our analysis shows that GW190814's source is a binary with an unequal mass ratio q = 0.112 +0.008 −0.009 , with individual source masses m 1 = 23.2 +1.1 −1.0 M and m 2 = 2.59 +0.08 −0.09 M , as shown in Figure 3 . A summary of the inferred source properties is given in Table 1 . We assume a standard flat ΛCDM cosmology with Hubble constant H 0 = 67.9 km s −1 Mpc −1 (Ade et al. 2016) .
We report detailed results obtained from the two precessing BBH signal models including subdominant multipole moments: Phenom PHM and EOBNR PHM. In order to compare the template models, we compute their Bayes factor (log 10 B). We find no significant evidence that one waveform family is preferred over the other as the Bayes factor between Phenom PHM and EOBNR PHM is log 10 B 1.0. As a result, we combine the posterior samples with equal weight, in effect marginalizing over a discrete set of signal models with a uniform probability. This is shown in the last column of Table 1 , and we refer to these values throughout the paper unless stated otherwise.
We find that the secondary mass lies in the range 2.50-2.67 M . This inferred secondary mass exceeds the bounds of the primary component in GW190425 (1.61-2.52 M ; Abbott et al. 2020a ) and the most massive known pulsar in the Galaxy: 2.14 +0.10 −0.09 M at 68.3% credible interval (Cromartie et al. 2019 ). Furthermore, the secondary is more massive than bounds on the maximum NS mass from studies of the remnant of GW170817, and from theoretical (Abbott et al. 1 In the co-precessing frame the EOBNR model includes the (l, m) = (2, ±2), (2, ±1), (3, ±3), (4, ±4) and (5, ±5) multipoles, and the Phenom model includes the (2, ±2), (2, ±1), (3, ±3) , (3, ±2) , (4, ±4) and (4, ±3) multipoles. (Farr & Chatziioannou 2020) . The inferred secondary mass is comparable to the putative BH remnant mass of GW170817 (Abbott et al. 2019b) .
The primary object is identified as a BH based on its measured mass of 23.2 +1.1 −1.0 M . Due to accurately observing the frequency evolution over a long inspiral, the chirp mass is well constrained to 6.09 +0.06 −0.06 M . The inferred mass ratio q = 0.112 +0.008 −0.009 makes GW190814 only the second gravitational-wave observation with a significantly unequal mass ratio (Abbott et al. 2019a (Abbott et al. , 2020d . Given that this system is in a region of the parameter space that has not been explored via gravitationalwave emission previously, we test possible waveform systematics by comparing the Phenom and EOB waveform families. Differences in the inferred secondary mass are shown in Figure 4 . The results indicate that the inferred secondary mass is robust to possible waveform systematics, with good agreement between the Phenom PHM and EOBNR PHM signal models. Signal models that exclude higher multipoles or precession do not constrain the secondary mass as well.
The time delay of a signal across a network of gravitational wave detectors, together with the relative amplitude and phase at each detector, allows us to measure the location of the GW source on the sky (Abbott et al. 2020b) . We localize GW190814's source to within 18.5 deg 2 at 90% probability, as shown in Figure 2 . This is comparable to the localization of GW170817 (Abbott et al. 2017a (Abbott et al. , 2019a .
Spins are a fundamental property of BHs. Their magnitude and orientation carry information regarding the evolution history of the binary. The effective inspiral spin parameter χ eff (Damour 2001; Racine 2008; Ajith et al. 2011; Santamaría et al. 2010) contains information about the spin components that are perpendicular to the orbital plane. We infer that χ eff = −0.002 +0.060 −0.061 . The tight constraints are consistent with being able to measure the phase evolution from the long inspiral.
Orbital precession occurs when there is a significant spin component in the orbital plane of the binary (Apostolatos et al. 1994) . We parameterize precession by the effective precession spin parameter 0 ≤ χ p ≤ 1 (Schmidt et al. 2015) . This effect is difficult to measure for face-on and face-off systems (Apostolatos et al. 1994; Buonanno et al. 2003; Vitale et al. 2014 Vitale et al. , 2017 Fairhurst et al. 2019a,b) . GW190814 constrains the inclination of the binary to be Θ = 0.8 +0.3 −0.2 rad. Since the system is neither face-on nor face-off, we are able to put strong constraints on the precession of the system: χ p = 0.04 +0.04 −0.03 . This is both the strongest constraint on the amount of precession for any gravitational-wave detection to date, and the first gravitational-wave measurement which conclusively measures near-zero precession (Abbott et al. 2019a (Abbott et al. , 2020a .
By computing the Bayes factor between a precessing and non-precessing signal model (log 10 B ∼ 0.5 in favor of precession), we find inconclusive evidence for in-plane spin. This is consistent with the inferred power from precession SNR ρ p (Fairhurst et al. 2019a,b) , whose recovered distribution resembles that expected in the absence of any precession in the signal; see Figure 5 . The ρ p calculation assumes a signal dominated by the = 2 mode; however, we have verified that the contribution of higher harmonics to the measurement of spin precession is subdominant by a factor of 5. The data are therefore consistent with the signal from a non-precessing system. Figure 4 shows that signal models including spinprecession effects give tighter constraints on the secondary mass compared to their non-precessing equivalents. Signal models that include spin-precession effects can constrain χ p , whereas non-precessing signal models cannot provide information on in-plane spin components. In all analyses, we assume a prior equivalent to spin orientations being isotropically distributed. We find that the data are inconsistent with large χ p and consistent with any secondary spin. Therefore, for precessing signal models the allowed q-χ eff parameter space is restricted, which helps to break the degeneracy (Poisson & Will 1995; Baird et al. 2013; Farr et al. 2016; Baird et al. 2013; Ng et al. 2018) . Consequently, the extra information from constraining χ p to small values enables a more precise measurement of the secondary mass.
The asymmetry in the masses of GW190814 means that the spin of the more massive object dominates contributions to χ eff and χ p . As both χ eff and χ p are tightly constrained, we are able to bound the primary spin of GW190814 to be χ 1 ≤ 0.07, as shown in Figure 6 . This is the strongest constraint on the primary spin for any gravitational-wave event to date (Abbott et al. 2019a (Abbott et al. , 2020a .
The joint posterior probability of the magnitude and orientation of χ 1 and χ 2 are shown in Figure 6 . Deviations from uniform shading indicate a spin property measurement. The primary spin is tightly constrained to small magnitudes, but its orientation is indistinguish- −0.02 , respectively. The final spin is lower than for previous mergers (Abbott et al. 2019a (Abbott et al. , 2020d , as expected from the low primary spin and smaller orbital contribution due to the asymmetric masses.
The relative importance of a subdominant multipole moment increases with mass ratio. Each subdominant multipole moment has a different angular dependence on the emission direction. With significant evidence for multipoles other than the dominant ( , m) = (2, 2) quadrupole, we gain an independent measurement of the inclination of the source. This allows for the distanceinclination degeneracy to be broken (Cutler & Flanagan 1994; Abbott et al. 2016f; Usman et al. 2019; Kalaghatgi et al. 2020) . Measuring higher-order multipoles therefore gives more precise measurements of source parameters (Van Den Broeck & Sengupta 2007a,b; Kidder 2008; Blanchet et al. 2008; Mishra et al. 2016; Kumar et al. 2019) .
GW190412 was the first event where there was significant evidence for higher-order multipoles (Payne et al. 2019; Kumar et al. 2019; Abbott et al. 2020d ). GW190814 exhibits stronger evidence for higher-order multipoles, with log 10 B 9.6 in favor of a higher- GW190814's stronger evidence for higher multipoles is expected given its more asymmetric masses and the larger network SNR.
The orthogonal optimal SNR of a subdominant multipole is calculated by decomposing each multipole into components parallel and perpendicular to the dominant harmonic (Mills & Fairhurst 2020; Abbott et al. 2020d) . We infer that the orthogonal optimal SNR of the ( , m) = (3, 3) multipole is 6.6 +1.3 −1.4 , as shown in Figure 5 . This is the strongest evidence for measuring a subdominant multipole to date (Payne et al. 2019; Kumar et al. 2019; Abbott et al. 2020d ).
Finally, we perform two complementary analyses involving time-frequency tracks in the data to provide further evidence for the presence of higher multipoles in the signal. In the first approach (also outlined in Abbott et al. 2020d , Section 4) we predict the time-frequency track of the dominant (2, 2) multipole in the LIGO Livingston detector (as seen in Figure 1 , middle panel) from an EOBNR HM parameter estimation analysis. This analysis collects energies along a time-frequency track which is α × f 22 (t), the (2, 2) multipole's instantaneous frequency, where α is a dimensionless parameter Abbott et al. 2020d) . We find prominent peaks in Y (α), the energy in the pixels along the α-th track defined in Abbott et al. (2020d) , at α = 1 and 1.5, as can be seen from the on-source curve in the top panel of Figure 7 . These peaks correspond to the m = 2 and m = 3 multipole predictions in the data containing the signal (on-source data). We also compute a detection statistic β (Roy et al. 2019) of 10.09 for the presence of the m = 3 multipole with a p-value of < 2.5 × 10 −4 , compared to a background distribution estimated over 18 hours of data adjacent to the event (off-source data), where the largest background β is 7.59. The significant difference between on-and off-source values provides much stronger evidence for the presence of higher multipoles than what is reported for GW190412 (Abbott et al. 2020d ).
The second analysis uses waveform-agnostic methods to reconstruct the signal. It then compares the observed coherent signal energy in the LIGO Hanford-LIGO Livingston-Virgo network of detectors, as identified by the cWB detection pipeline (Klimenko et al. 2016) , with the predictions of a waveform model without higher multipoles (EOBNR; Prodi et al. 2020) to investigate if the description of the underlying signal is incomplete if we do not include contributions from the m = 3 multipole in our waveform model. We compute a test statistic, the squared sum of the coherent residuals estimated over selected time-frequency tracks parameterized in terms of the same α parameter defined in the previous analysis Abbott et al. 2020d) . Each time-frequency track centered on α includes frequencies within [α − 0.1, α + 0.1] × f 22 (t), and times within [t merger − 0.5 s, t merger − 0.03 s], where f 22 (t) and t merger correspond to the maximum likelihood template from the EOBNR parameter estimation analysis. We further compute a background distribution using simulated signals in off-source data (Prodi et al. 2020) , and compute p-values for the on-source results as a function of α (Figure 7, bottom panel) . We find a minimum p-value of 6.8 × 10 −3 at α = 1.5, providing strong evidence that the disagreement between the actual event and the EOBNR prediction is because of the absence of the m = 3 multipoles in the waveform model. The local minimum near α = 2 is not an indication of the m = 4 multipoles, but rather a statistical fluctuation which is consistent with similar behaviour seen for studies with simulated signals described in detail in Prodi et al. (2020) .
Although the two time-frequency analyses are similar in motivation, the latter differs from the former in that it is not restricted to data from just one detector, but rather uses the coherent signal energy across the three-detector network. Both analyses point to strong evidence for the presence of higher multipoles in the signal.
GW190814 is the gravitational-wave event with the most unequal mass ratio to date, and can therefore be used to test general relativity (GR) in a region of parameter space previously unexplored with strong-field tests of GR (Abbott et al. 2016g, 2019e,d) . The asymmetric nature of a system excites the higher multipole moments of the gravitational signal, which allows us to test the multipolar structure of gravity (Kastha et al. 2018 (Kastha et al. , 2019 Dhanpal et al. 2019; Islam et al. 2020 ). The addition of information from the higher harmonics of a signal also breaks certain degeneracies in the description of the source, and could potentially enable us to place stronger constraints on certain deviations from GR (Van Den Broeck & Sengupta 2007b,a). We perform several null tests of GR using GW190814. These tests assume GW190814 is a (quasi-circular) BBH merger as described in GR, and look for inconsistencies between the observed signal and predictions of the theory. An inconsistency might arise from an incomplete understanding of the underlying signal (or noise), and could indicate a non-BBH nature of the signal or a potential departure from GR.
First, as a consistency test of the signal reconstruction, we subtract from the data the maximum likelihood compact binary coalescence waveforms, Phenom (Khan et al. 2016) , Phenom HM (Kalaghatgi et al. 2020) , Phenom PHM , and EOBNR PHM (Ossokine et al. 2020 ) and analyze 4 s of the resulting residual data centered around the time of merger with the morphology-independent transient analysis BayesWave . We measure the 90% credible upper limit on the coherent SNR, ρ 90 , and compare it to the SNR, ρ N 90 , recovered by analyzing 175 randomly selected data segments in surrounding time (off-source data) with the same configuration settings. If the residual data are consistent with the noise, we expect ρ 90 to be consistent with ρ N 90 . We compute the p-value by comparing the distribution of ρ N 90 to ρ 90 through p = P (ρ N 90 < ρ 90 ). We obtain p-values of 0.59, 0.82, 0.82, and 0.75 for Phenom, Phenom HM, Phenom PHM, and EOBNR PHM, respectively. Hence, we find no evidence for deviations in the behavior of the residual data stream.
We also look for deviations in the spin-induced quadrupole moments of the binary components. According to the no-hair conjecture (Carter 1971; Hansen 1974) the multipole moments of a Kerr BH are completely described by its mass and spin angular momentum. At leading order in spin, the spin-induced quadrupole moment scalar is (Hartle 1967; Pappas & Apostolatos 2012) , Q = −κa 2 m 3 , where (m, a) are the mass and dimensionless spin of the compact object, and κ is a dimensionless deformation parameter characterizing deviations in the spin-induced quadrupole moment. Kerr BHs have κ = 1 (Thorne 1980) , while κ ∼ 2-14 for NSs (depending on the equation of state) and κ ∼ 10-150 for spinning boson stars with large self-interaction (Ryan 1997) . The deformation parameter can even be negative for (slowly-rotating, thin-shelled) gravastars (Uchikata et al. 2016) . Hence, an accurate measurement of κ sheds light onto the nature of the compact object. For compact binaries, the spin-induced quadrupole moment terms appear at second post-Newtonian order (Poisson 1998) . For Kerr BHs in GR, κ 1 = κ 2 = 1, where κ 1 , κ 2 are the individual deformation parameters of the primary and secondary compact objects in the binary. Since κ 1 and κ 2 are strongly degenerate in the gravitational waveform, we instead measure a linear symmetric combination of these quantities, κ s = (κ 1 +κ 2 )/2, which is 1 for a BBH in GR. The posteriors on κ s are relatively uninformative, and nearly span the prior range of [0, 500], with increased support at κ s = 0 relative to the prior. The upper bound of the prior was chosen to accommodate all the objects listed above. The result shows that GW190814 is consistent with having a BBH source described by GR. However, the broad posterior means that we cannot exclude the possibility that one or both components of the source is not a BH. We can attempt to understand this result in terms of the spin measurements for the binary. The measurements of κ s and a non-zero χ eff are highly correlated (Krishnendu et al. 2019) , and for a system with small χ eff the bounds on the measured value of κ s are weak.
Finally, we investigate the source dynamics of the binary through a parameterized test of gravitational waveform generation, where we allow for the coefficients describing the post-Newtonian inspiral of a BBH coalescence to deviate away from their predictions in GR (Arun et al. 2006a,b; Yunes & Pretorius 2009; Mishra et al. 2010; Cornish et al. 2011; Li et al. 2012; Meidam et al. 2018) . We use an aligned-spin EOB waveform without higher modes (EOBNR), and find no deviations in the post-Newtonian coefficients from their nominal values in GR. In summary, none of our tests of GR indicate any departure from the predictions of the theory, and GW190814 is consistent with the description of a compact binary merger in GR.
The highly unequal mass ratio of 0.112 +0.008 −0.009 and unusual secondary mass of 2.59 +0.08 −0.09 M make the source of GW190814 unlike any other compact binary coalescence observed so far. The average mass ratio for BBH coalescences detected by the LVC during O1 and O2 is 0.9 (Roulet & Zaldarriaga 2019) , and an inference of the underlying population predicted that 99% of detectable BBHs have mass ratios q ≥ 0.5 (Fishbach & Holz 2020 ). However, the paucity of events from O1 and O2 means that this picture is limited. Indeed, the discovery of GW190412 has already changed the picture substantially (Abbott et al. 2020d ).
GW190814's secondary mass lies in the hypothesized lower mass gap of 2. 5-5 M (Bailyn et al. 1998; Özel et al. 2010; Farr et al. 2011; Özel et al. 2012 ) between known NSs and BHs. It is heavier than the most massive pulsar in the Galaxy (Cromartie et al. 2019) , and almost certainly exceeds the mass of the 1.61-2.52 M primary component of GW190425, which is itself an outlier relative to the Galactic population of BNSs (Abbott et al. 2020a ). On the other hand, it is comparable in mass to two BH candidates: the 2.7 M merger remnant of GW170817 (Abbott et al. 2019b ) and the 2.6-6.1 M compact object (95% confidence interval) discovered by Thompson et al. (2019) . 2 It is also comparable to the millisecond pulsar PSR J1748−2021B (Freire et al. 2008) , whose mass is claimed as 2.74 +0.21 −0.21 M at 68% confidence. However, this estimate, obtained via measurement of the periastron advance, could be inaccurate if the system inclination is low or the pulsar's companion is rapidly rotating (Freire et al. 2008) . In sum, it is not clear if GW190814's secondary is a BH or a NS.
GW190814 poses a challenge for our understanding of the population of merging compact binaries. In what follows, we estimate the merger rate density of the compact binary subpopulation represented by this source, investigate the nature of its secondary component and possible implications for the NS equation of state, discuss how the system may have formed, and study its implications for cosmology.
Given the unprecedented combination of component masses found in GW190814, we take the system to represent a new class of compact binary mergers, and use our analysis of its source properties to estimate a merger rate density for GW190814-like events. Following a method described in Kim et al. (2003) , we calculate a simple, single-event rate density estimate R according to our sensitivity to a population of systems drawn from the parameter-estimation posteriors. As in Abbott et al. (2020a) , we calculate our surveyed spacetime volume V T semi-analytically, imposing singledetector and network SNR thresholds of 5 and 10, respectively (Tiwari 2018) . The semi-analytic V T for GW190814 is then multiplied by a calibration factor to match results from the search pipelines assuming a onceper-century FAR threshold. The sensitivity of a search pipeline is estimated using a set of simulated signals. For computational efficiency, this was done using preexisting search pipeline simulations and the mass properties were not highly optimized. However, given that we are estimating a rate based on a single source, the cal-ibration errors are much smaller than the statistical errors associated with the estimate. The simulated sources were uniformly distributed in comoving volume, component masses, and component spins aligned with the orbital angular momentum. For O1 and O2, the simulated BH mass range was 5-100 M , but for the first part of O3 we are analyzing here, the injected range was 2.5-40 M (following our updated knowledge of the BH mass distribution); the NS mass range was 1-3 M , and component spins are < 0.95. As GW190814 occurred when LIGO Hanford was not in nominal observing mode, it is not included in the production PyCBC results, and we use GstLAL results to calculate the merger rate.
We assume a Poisson likelihood over the astrophysical rate with a single count and we apply a Jeffreys R −1/2 prior to obtain rate posteriors. The analysis was done using samples from the Phenom PHM posterior and separately from the EOBNR PHM posterior, producing the same result in both cases. We find the merger rate density of GW190814-like systems to be 7 +16 −6 Gpc −3 yr −1 . As a consistency check, we used the PyCBC search results to calculate an upper limit. Repeating the rate calculation with a PyCBC-based V T calibration and zero event count, we obtain an upper limit consistent (to within 10%) with the upper limit of the merger rate estimated using GstLAL search results. We conclude that the uncertainty in our estimate of the rate density for the class of mergers represented by GW190814 is primarily dominated by Poisson statistics.
The primary mass measurement of 23.2 +1.1 −1.0 M securely identifies the heavier component of GW190814 as a BH, but the secondary mass of 2.59 +0.08 −0.09 M may be compatible with either a NS or a BH depending on the maximum mass supported by the unknown NS equation of state (EOS). The source's asymmetric masses, the non-detection of an electromagnetic counterpart and the lack of a clear signature of tides or spin-induced quadrupole effects in the waveform do not allow us to distinguish between a BBH or a NSBH. Instead, we rely on comparisons between m 2 and different estimates of the maximum NS mass, M max , to indicate the source classification preferred by data: if m 2 > M max , then the NSBH scenario is untenable.
While some candidate EOSs from nuclear theory can support nonrotating NSs with masses of up to ∼ 3 M (e.g., Müller & Serot 1996) , such large values of M max are disfavored by the relatively small tidal deformabilities measured in GW170817 (Abbott et al. 2017a (Abbott et al. , 2019b , which correlate with smaller internal pressure gradients as a function of density and hence a lower threshold for gravitational collapse. By adopting a phenomenological model for the EOS, conditioning it on GW170817, and extrapolating the constraints to the high densities relevant for the maximum mass, Lim & Holt (2019) and Essick et al. (2020) place M max 2.3 M . Similarly, the EOS inference reported in Abbott et al. (2018) , based on an analysis of GW170817 with a spectral parameterization (Lindblom 2010; Lindblom & Indik 2012 for the EOS, implies a 90% credible upper bound of M max ≤ 2.43 M , with tenuous but non-zero posterior support beyond 2.6 M . We calculate the corresponding M max posterior distribution, shown in the right panel of Figure 3 , from the GW170817informed spectral EOS samples used in Abbott et al. (2018) by reconstructing each EOS from its parameters and computing its maximum mass. Comparison with the m 2 posterior suggests that the secondary component of GW190814 is probably more massive than this prediction for M max : the posterior probability of m 2 ≤ M max , marginalized over the uncertainty in m 2 and M max , is only 3%. Nevertheless, the maximum mass predictions from these kinds of EOS inferences come with important caveats: their extrapolations are sensitive to the phenomenological model assumed for the EOS; they use hard M max thresholds on the EOS prior to account for the existence of the heaviest Galactic pulsars, which is known to bias the inferred maximum mass distribution towards the threshold (Miller et al. 2020) ; and they predate the NICER observatory's recent simultaneous mass and radius measurement for J0030+0451, which may increase the M max estimates by a few percent ) because it favors slightly stiffer EOSs than GW170817 (Raaijmakers et al. 2019; Riley et al. 2019; Miller et al. 2019; Jiang et al. 2020) .
NS mass measurements also inform bounds on M max independently of EOS assumptions. Fitting the known population of NSs in binaries to a double-Gaussian mass distribution with a high-mass cutoff, Alsing et al. (2018) obtained an empirical constraint of M max ≤ 2.6 M (one-sided 90% confidence interval). Farr & Chatziioannou (2020) recently updated this analysis to include PSR J0740+6620 (Cromartie et al. 2019 ), which had not been discovered at the time of the original study. Based on samples from the Farr & Chatziioannou (2020) maximum-mass posterior distribution, which is plotted in the right panel of Figure 3 , we find M max = 2.25 +0.81 −0.26 M . In this case, the posterior probability of m 2 ≤ M max is 29%, again favoring the m 2 > M max scenario, albeit less strongly because of the distribution's long tail up to ∼ 3 M . However, the empirical M max prediction is sensitive to selection effects that could potentially bias it (Alsing et al. 2018 ). In particular, masses are only measurable for binary pulsars, and the mass distribution of isolated NSs could be different. Additionally, the discovery of GW190425 (Abbott et al. 2020a) should also be taken into account in the population when predicting M max .
Finally, the NS maximum mass is constrained by studies of the merger remnant of GW170817. Although no postmerger gravitational waves were observed (Abbott et al. 2017g, 2019f) , modeling of the associated kilonova (Abbott et al. 2017b; Kasen et al. 2017; Villar et al. 2017; Cowperthwaite et al. 2017; Abbott et al. 2017d) suggests that the merger remnant collapsed to a BH after a brief supramassive or hypermassive NS phase during which it was stabilized by uniform or differential rotation. Assuming this ultimate fate for the merger remnant immediately implies that no NS can be stable above ∼ 2.7 M , but it places a more stringent constraint on NSs that are not rotationally supported. The precise mapping from the collapse threshold mass of the remnant to M max depends on the EOS, but by developing approximate prescriptions based on sequences of rapidly rotating stars for a range of candidate EOSs, M max has been bounded below approximately 2.2-2.3 M (Margalit & Metzger 2017; Rezzolla et al. 2018; Ruiz et al. 2018; Shibata et al. 2019; Abbott et al. 2020c) . Although the degree of EOS uncertainty in these results is difficult to quantify precisely, if we take the more conservative 2.3 M bound at face value, then m 2 is almost certainly not a NS: the m 2 posterior distribution has negligible support below 2.3 M .
Overall, these considerations suggest that GW190814 is probably not the product of a NSBH coalescence, despite its preliminary classification as such. Nonetheless, the possibility that the secondary component is a NS cannot be completely discounted due to the current uncertainty in M max .
There are two further caveats to this assessment. First, because the secondary's spin is unconstrained, it could conceivably be rotating rapidly enough for m 2 to exceed M max without triggering gravitational collapse: rapid uniform rotation can stabilize a star up to ∼ 20% more massive than the nonrotating maximum mass (Cook et al. 1994) , in which case only the absolute upper bound of ∼ 2.7 M is relevant. However, it is very unlikely that a NSBH system could merge before dissipating such extreme natal NS spin angular momentum.
Second, our discussion has thus far neglected the possibility that the secondary component is an exotic compact object, such as a boson star (Kaup 1968) or a gravastar (Mazur & Mottola 2004) , instead of a NS or a BH. Depending on the model, some exotic compact objects can potentially support masses up to and beyond Figure 8 . Constraints on the NS EOS assuming GW190814 was produced by a BBH (blue) or a NSBH (orange) coalescence. The 90% and 50% credible contours of the posterior in the pressure-density plane are shown. The constraints are calculated by assuming a spectral decomposition for the EOS, following Abbott et al. (2018) . The BBH constraints are identical to those from the analysis of GW170817, while for the NSBH case the posterior is reweighted by the probability that each EOS's maximum mass is at least m2. The dashed lines indicate the 90% credible region of the prior.
Since the NSBH scenario cannot be definitively ruled out, we examine GW190814's potential implications for the NS EOS, assuming that the secondary proves to be a NS. This would require M max to be no less than m 2 , a condition that severely constrains the distribution of EOSs compatible with existing astrophysical data. The combined constraints on the EOS from GW170817 and this hypothetical maximum mass information are shown in Figure 8 . Specifically, we have taken the spectral EOS distribution conditioned on GW170817 from Abbott et al. (2018) and reweighted each EOS by the probability that its maximum mass is at least as large as m 2 . The updated posterior favors stiffer EOSs, which translates to larger radii for NSs of a given mass. The corresponding constraints on the radius and tidal deformability of a canonical 1.4 M NS are R 1.4 = 12.9 +0.8 −0.7 km and Λ 1.4 = 616 +273 −158 .
The source of GW190814 represents a previously undetected class of coalescences that has the potential to shed light on the formation of merging compact-object binaries with highly asymmetric masses.
Electromagnetic observations of Galactic NSs and stellar-mass BHs suggest a dearth of compact objects in the ∼ 2.5 M to 5 M range (Bailyn et al. 1998; Özel et al. 2010; Farr et al. 2011; Özel et al. 2012) . Observa-tions of a few candidates with masses in this range seem to disfavor the existence of a gap (Freire et al. 2008; Neustroev et al. 2014; Giesers et al. 2018; Wyrzykowski & Mandel 2020; Thompson et al. 2019 ), but whether the mass gap is physical or caused by selection biases is still a matter of debate (e.g., Kreidberg et al. 2012) .
From a theoretical point of view, accurately calculating the masses of compact remnants at formation is challenging, because it depends on the complex physics of the supernova explosion and the details of stellar evolution, especially for the late evolutionary stages of massive stars (Janka 2012; Müller 2016; Burrows et al. 2018; Burrows et al. 2019) . Whether the models favor the presence of a gap or a smooth transition between NSs and BHs is still unclear, and in fact some models have been developed with the purpose of reproducing this lower mass gap (Ugliano et al. 2012; Fryer et al. 2012; Kochanek 2014; Sukhbold & Woosley 2014; Ertl et al. 2016) . Therefore, our robust discovery of an object with a well-constrained mass in this regime may provide crucial constraints on compact-object formation models. In fact, GW190814 demonstrates the need to adjust remnant mass prescriptions previously designed to produce a perceived mass gap. The combination of mass ratio and component masses challenges most results obtained from population synthesis simulations for isolated binaries (Dominik et al. 2012 (Dominik et al. , 2015 Marchant et al. 2017; Kruckow et al. 2018; Neijssel et al. 2019; Mapelli et al. 2019; Spera et al. 2019; Olejak et al. 2020) .
Population synthesis models distinguish between NSs and BHs using only a mass threshold, which is generally in the range 2-3 M . Thus, depending on the adopted threshold and on the adopted supernova explosion model, a GW190814-like event may be labeled as either a NSBH merger or a BBH merger. Most BBH mergers have q > 0.5, while the distributions of merging NSBH binaries suggest that systems with q 0.1 may be up to ∼ 10 3 times less common than more symmetric ones (q > 0.1) and that the mass-ratio distribution peaks at q ≈ 0.2. Furthermore, models tend to favor mergers of massive ( 1.3 M ) NSs with relatively small BHs ( 15 M ) in environments with sub-solar metallicity (Z 0.5 Z ). The tendency to disfavor mergers with highly asymmetric masses in isolated binaries may be the consequence of mass transfer (e.g., Postnov & Yungelson 2014) and common envelope episodes (e.g., Ivanova et al. 2013 ) that cause systems with initially asymmetric masses to evolve towards more symmetric configurations. Overall, producing mergers with such unequal masses, with a secondary in the perceived mass gap, and at the rate implied by this discovery is a challenge for current models.
Nevertheless, particular choices of poorly constrained assumptions within rapid population synthesis models may increase the number of mergers with q 0.1 so that the latter may be only a few times less common than (or even comparable to) systems with q 0.2 (e.g., Eldridge & Stanway 2016; Eldridge et al. 2017; .
Another possibility is that GW190814 is of dynamical origin. Dynamical exchanges in dense stellar environments tend to pair up massive compact objects with similar masses (e.g., Sigurdsson & Hernquist 1993) . This process is effective for globular clusters, where compactobject binaries may undergo tens of exchanges before they get ejected from the cluster (Portegies Zwart & McMillan 2000; Rodriguez et al. 2016; Park et al. 2017; Askar et al. 2017; Rodriguez et al. 2019 ). For such environments, models predict that most merging BBHs have q 1 (e.g., Rodriguez et al. 2016) , and the formation of NSBH binaries is highly suppressed because BHs dynamically dominate the cores over the complete lifetime of the clusters, preventing the interactions between BHs and NSs, with the consequence that the merger rate of NSBH binaries in globular clusters in the local Universe is ∼ 10 −2 -10 −1 Gpc −3 yr −1 (Clausen et al. 2013; Ye et al. 2020; Arca Sedda 2020) . The rate for GW190814-like events, with a secondary in the perceived mass gap, is likely even lower. In contrast, the NSBH merger rate may be significantly higher in young star clusters (e.g., Ziosi et al. 2014 ) and the latter can effectively increase the number of progenitors leading to merging compactobject binaries with q 0.15 (Di Carlo et al. 2019; Rastello et al. 2020) . Thus, young star clusters may be promising hosts for GW190814-like events, but the parameter space relevant for GW190814 is mostly unexplored in the context of star clusters.
In dense stellar environments, GW190814-like systems may also form from a low-mass merger remnant that acquires a BH companion via dynamical interactions (Gupta et al. 2020) . Gupta et al. (2020) predicts a population of second-generation BHs in the 2.2-3.8 M range, with a peak in the distribution at 2.6 M , assuming a double-Gaussian mass distribution for the NSs. However, recent dynamical simulations of globular clusters (e.g., Ye et al. 2020) find the subsequent merger of such a second-generation BH with a larger stellarmass BH to be exceedingly rare. A high component spin could be a distinguishing feature of a second-generation compact object, but the uninformative spin posterior for the lighter component of GW190814 provides no evidence for or against this hypothesis.
A GW190814-like merger may also have originated from a hierarchical triple in the field (e.g., Silsbee & Tremaine 2017; Fragione & Loeb 2019; Antonini et al. 2017) , from a wide hierarchical quadruple system (Safarzadeh et al. 2020) , or from hierarchical triples in galactic centers, where the tertiary body is a supermassive BH (Antonini & Perets 2012; Stephan et al. 2019; Petrovich & Antonini 2017; Hoang et al. 2018; . Specifically, Safarzadeh et al. (2020) explore the possibility that a second-generation remnant with mass 3 M may merge with a 30 M BH, catalyzed by a 50 M -BH perturber. The mass-ratio distributions of BBH and NSBH mergers from hierarchical systems are similar to those of field binaries and it is unclear whether hierarchies may enhance the formation of merging compact-object binaries with highly asymmetric masses (e.g., Silsbee & Tremaine 2017) .
Disks of gas around supermassive BHs in active galactic nuclei may be promising environments for the formation of GW190814-like systems. For such environments, theoretical models show that merging compact-object binaries with asymmetric masses are likely, but cannot necessarily accommodate masses as low as the secondary mass of GW190814 (e.g., Yang et al. 2019 ). However, McKernan et al. (2020) show that the median mass ratio of NSBH mergers in active galactic nucleus disks may be as low as ∼ 0.07.
We conclude that the combination of masses, mass ratio and inferred rate of GW190814 is challenging to explain, but potentially consistent with multiple formation scenarios. However, it is not possible to assess the validity of models that produce the right properties but do not make quantitative predictions about formation rates, even at some order-of-magnitude level.
Young star clusters and active galactic nucleus disks seem to be more promising hosts for GW190814-like mergers, since both these environments may enhance the formation of either progenitors of or directly merging compact-object binaries with more asymmetric masses to relevant rates. In contrast, globular-cluster models provide more robust predictions, showing that GW190814-like mergers with such asymmetric masses are outliers in the population predictions, even though a revision of the remnant-mass prescription is still needed. Isolated binaries binaries could prove possible progenitors provided similar revisions are implemented. The importance of field multiples remains to be fully explored. Future gravitational-wave observations will provide further insights into the dominance of different channels.
Luminosity distances inferred directly from observed gravitational-wave events can be used with measurements of source redshifts in the electromagnetic spectrum to constrain cosmological parameters (Schutz 1986 ). Redshifts can be either obtained directly from counterparts to the gravitational-wave source (Holz & Hughes 2005) , as was the case for GW170817 (Abbott et al. 2017a,b,e) , by cross-correlation of the gravitational-wave localization posterior with catalogs of galaxy redshifts (Del Pozzo 2012; Nair et al. 2018; Chen et al. 2018; Fishbach et al. 2019; Gray et al. 2020; Soares-Santos et al. 2019; Abbott et al. 2019c ), by exploiting information in the neutron star equation of state (Messenger & Read 2012) , or by using the redshifted masses inferred from the gravitational wave observation and assumptions about the mass distribution of the sources (Chernoff & Finn 1993; Farr et al. 2019) . At current sensitivities, the cosmological parameter to which LIGO-Virgo observations are most sensitive is the Hubble constant, H 0 . The gravitational-wave observation of GW170817 provided a posterior on H 0 with mode and 68.3% highest posterior density interval of H 0 = 69 +22 −8 km s −1 Mpc −1 (Abbott et al. 2017e, 2019b , assuming a flat prior on H 0 .
GW190814 is the best localized dark siren, i.e., gravitational-wave source without an electromagnetic counterpart, observed to date, and so it is a good candidate for the statistical cross-correlation method. For a fixed reference cosmology (Ade et al. 2016) , the GLADE galaxy catalog (Dálya et al. 2018 ) is approximately 40% complete at the distance of GW190814 and contains 472 galaxies within the 90% posterior credible volume of GW190814. To obtain a constrain on H 0 , we use the methodology described in Abbott et al. (2019c) and the GLADE catalog. We take a flat prior for H 0 ∈ [20, 140] km s −1 Mpc −1 and assign a probability to each galaxy that it is the true host of the event that is proportional to its B-band luminosity. Using the posterior distribution on the distance obtained from the combined PHM samples, we obtain H 0 = 75 +59 −13 km s −1 Mpc −1 using GW190814 alone (mode and 68.3% highest posterior density interval; the median and 90% symmetric credible interval is H 0 = 83 +55 −53 km s −1 Mpc −1 ), which can be compared to H 0 = 75 +40 −32 km s −1 Mpc −1 (Soares-Santos et al. 2019) obtained using the dark siren GW170814 alone. The GW190814 result is the most precise measurement from a single dark siren observation to date, albeit comparable to the GW170814 result, which is expected given GW190814's small localization volume (∼ 39000 Mpc 3 ). The result is not very constraining, with the 68.3% highest posterior density interval com-prising 60% of the prior range. Combining the result for GW190814 with the result obtained from GW170817, we see an improvement over the GW170817-only result, to H 0 = 70 +17 −8 km s −1 Mpc −1 (the median and 90% symmetric credible interval is H 0 = 77 +33 −23 km s −1 Mpc −1 ). This result is not yet sufficiently constraining to provide further insight into current tensions in low and high redshift measurements of the Hubble constant (Verde et al. 2019 ), but these constraints will continue to improve as further gravitational-wave observations are included (e.g., projections in Chen et al. 2018; Vitale & Chen 2018; Gray et al. 2020; Feeney et al. 2019) .
During their third observing run, on 2019 August 14, 21:10:39 UTC, the LIGO and Virgo detectors observed GW190814, a novel source unlike any other known compact binary coalescence. Thanks in part to the observation of significant power in subdominant multipoles of the gravitational radiation, and the conclusive measurement of little to no spin precession, we obtain precise measurements of its physical source properties that clearly set it apart from other compact binaries.
In particular, (a) its mass ratio of q = 0.112 +0.008 −0.009 is the most unequal ever observed with gravitational waves, (b) the bound χ 1 ≤ 0.07 on the spin of the 23.2 +1.1 −1.0 M BH is the strongest constraint on a primary spin for any gravitational-wave source to date, and (c) the secondary mass measurement of 2.59 +0.08 −0.09 M makes it the lightest BH or the heaviest NS discovered in a double compact-object system. We find no evidence of measurable tidal effects in the signal, and no electromagnetic counterpart to the gravitational waves has been identified.
Comparisons between the secondary mass and several current estimates of the maximum NS mass suggest that GW190814 is unlikely to originate in a NSBH coalescence. Nevertheless, the M max estimates are uncertain enough that improved knowledge of the NS EOS or further observations of the astrophysical population of compact objects could alter this assessment. For this reason, we cannot firmly exclude the possibility that the secondary is a NS, nor can we be certain that it is a BH. Regardless, this event sheds new light on the compactobject mass distribution at the interface between known NSs and BHs.
The unique combination of masses and inferred merger rate for this event is difficult to produce in population synthesis models of multi-component systems in galactic fields or dense stellar environments. The discovery of GW190814 may therefore reshape our understanding of the processes by which the lightest BHs or the most massive NSs form. Based on our rate density estimate, we may reasonably expect to detect more systems of this kind after a year at design sensitivity. This discovery may prove to be the first hint of a larger population that could change our perspective on the formation and mass spectrum of compact objects.
Segments of data containing the signal from all three interferometers, and samples from the posterior distributions for the source parameters, are available from the Gravitational Wave Open Science Center (https://doi.org/10.7935/zzw5-ak90).
The software packages used in our analysis are open source.
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 Council of Scientific and 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.
Software: The detection of the signal and subsequent significance evaluation were performed with the GstLAL-based inspiral software pipeline (Cannon et al. 2012; Privitera et al. 2014; Messick et al. 2017; Sachdev et al. 2019; Hanna et al. 2020) , built on the LALSuite software library (LIGO Scientific Collaboration 2018), and with the PyCBC (Nitz et al. 2018 (Nitz et al. , 2019 Usman et al. 2016) and MBTAOnline (Adams et al. 2016) packages. Parameter estimation was performed with the LALInference (Veitch et al. 2015) and LALSimulation libraries within LALSuite (LIGO Scientific Collaboration 2018), as well as the Bilby and pBilby Libraries Smith & Ashton 2019 ) and the dynesty nested sampling package (Speagle 2020) . Interpretation and curation of the posterior samples was handled by the PESummary library (Hoy & Raymond 2020) . Estimates of the noise spectra were obtained using BayesWave . Plots were prepared with Matplotlib (Hunter 2007) . The sky map plot also used Astropy (http://www.astropy.org) a community-developed core Python package for Astronomy (Astropy Collaboration et al. 2013; Price-Whelan et al. 2018 ) and ligo.skymap (https://lscsoft.docs.ligo.org/ligo.skymap).