GW190412: Observation of a Binary-Black-Hole Coalescence with Asymmetric Masses

We report the observation of gravitational waves from a binary-black-hole coalescence during the first two weeks of LIGO's and Virgo's third observing run. The signal was recorded on April 12, 2019 at 05:30:44 UTC with a network signal-to-noise ratio of 19. The binary is different from observations during the first two observing runs most notably due to its asymmetric masses: a ~30 solar mass black hole merged with a ~8 solar mass black hole companion. The more massive black hole rotated with a dimensionless spin magnitude between 0.17 and 0.59 (90% probability). Asymmetric systems are predicted to emit gravitational waves with stronger contributions from higher multipoles, and indeed we find strong evidence for gravitational radiation beyond the leading quadrupolar order in the observed signal. A suite of tests performed on GW190412 indicates consistency with Einstein's general theory of relativity. While the mass ratio of this system differs from all previous detections, we show that it is consistent with the population model of stellar binary black holes inferred from the first two observing runs.

To be able to characterize the full range of potential systems, models of the gravitational radiation emitted by BBHs are continuously being improved. In particular, physical effects associated with unequal masses and misaligned spins have recently been extended in models covering the inspiral, merger and ringdown of BBHs [21][22][23][24][25][26][27][28][29][30][31][32][33]. For GW signals with sufficient signal-to-noise ratio (SNR), the inclusion of these effects is important to accurately infer the source parameters. In addition, improved signal models allow for stronger tests of the validity of general relativity (GR) as the correct underlying theory.
In this paper we report the detection of GWs from a BBH whose properties make it distinct from all other BBHs reported previously from the first two observing runs. The event, called GW190412, was observed on April 12, 2019 at 05:30:44 UTC by the Advanced Virgo detector and both Advanced LIGO detectors. While the inferred individual masses of the coalescing black holes (BHs) are each within the range of masses that have been observed before [7,[9][10][11][12], previously detected binaries all had mass ratios q = m 2 /m 1 (with m 1 ≥ m 2 ) that were consistent with unity [34]. GW190412, however, is the first observation of GWs from a coalescing binary with unequivocally unequal masses, q = 0.28 +0.13 −0.06 (median and 90% symmetric credible interval). The mass asymmetry of the system provides a second novelty of GW190412: the GWs carry subtle, but for the first time clearly measurable, imprints of radiation that oscillates not only at the binary's dominant GW emission frequency, but also at other frequencies with subdominant contributions. We introduce the nature of these corrections and present the source parameters inferred for GW190412 using signal models that include higher multipoles.
This paper is organized as follows: in Sec. II we report details on the detection of GW190412. The source properties are discussed in Sec. III. Sec. IV presents a suite of analyses illustrating that the observed data indeed contain measurable imprints of higher multipoles. In Sec. V we present tests of GR performed in this previously unexplored region of the parameter space. Implications for our understanding of the BBH population and formation channels are discussed in Secs. VI and VII.

II. DETECTORS & DETECTION
The third observing run of LIGO [35] and Virgo [14] began on 1 April 2019, and GW190412 occurred in the second week of the run. At the time of the event, both LIGO detectors and the Virgo detector were online and operating stably for over 3.5 hours. Strain data from around the time of GW190412 for all three detectors is shown in Figure 1, with excess power consistent with the observed signal present in all detectors. The relative sensitivity of the LIGO and Virgo detectors accounts for the difference in strength of the signal in the data.
LIGO and Virgo interferometers are calibrated using electrostatic fields and radiation pressure from auxiliary lasers at a known frequency and amplitude [36][37][38]. At the time of the event, the maximum calibration error at both LIGO sites is 7.0% in amplitude and 3.8 • in phase. At Virgo, the errors are 5.0% in amplitude and 7.5 • in phase.
Numerous noise sources that limit detector sensitivity are measured and subtracted, including noise from calibration lines and noise from the harmonics of the power mains. Similar to procedures from the second observing run [39,40], these noise sources are linearly subtracted from the data using auxiliary witness sensors. In O3, this procedure was completed as a part of the calibration pipeline [37], both in low latency and offline. Additional noise contributions due to non-linear coupling of the 60 Hz power mains are subtracted for offline analyses using coupling functions that rely on machine learning techniques [41]. GW190412 was initially detected by the GstLAL [43], SPIIR [44], cWB [45], MBTA [46], and PyCBC Live [47] pipelines running in low-latency configuration, and reported under the identifier S190412m. The GstLAL, SPIIR, MBTA, and PyCBC Live pipelines identify GW signals by matched-filtering [48][49][50] data using a bank of filter waveforms that cover a wide range of source parameters [51][52][53][54][55][56][57]. The coherent, unmodelled cWB pipeline [45], identifies clusters of coherent excess power with increasing frequency over time in data processed with the wavelet transform [58].
All analysis pipelines running in low-latency identified GW190412 as a confident event. The observed SNR from the GstLAL pipeline was 8.6 in LIGO Hanford, 15.6 in LIGO Livingston, and 3.7 in Virgo. GW190412 was identified with consistent SNR across all low-latency pipelines.
A GCN alert [59] announcing the event was publicly distributed 60 minutes after GW190412 and included an initial sky localization computed using a rapid Bayesian algorithm, BAYESTAR [60], applied to data from all available detectors. This sky localization constrained the position of the event to a 90% credible area of 156 deg 2 .
Additional offline analysis of the data from April 8 to April 18 was completed by the GstLAL [61,62] and Py-CBC [63,64] matched filtering pipelines, and the coherent, template independent cWB pipeline [45]. The offline analyses utilize an updated version of the calibration of the LIGO data [37] and additional data quality vetoes [65]. The GstLAL pipeline incorporates data from all three detectors, while the PyCBC and cWB pipelines only use data from the two LIGO detectors.
All three offline pipelines identified GW190412 as a highly significant GW event. This event was assigned a false-alarm rate (FAR) of < 1 per 1 × 10 5 years by the GstLAL pipeline and < 1 per 3 × 10 4 years by the Py-CBC pipeline. The template independent cWB pipeline assigned this event a FAR of < 1 per 1 × 10 3 years. As GW190412 was identified as more significant than any event in the computed background, the FARs assigned by all offline pipelines are upper bounds.
Validation procedures similar to those used to evaluate previous events [7,66] were used in the case of GW190412 to verify that instrumental artifacts do not affect the analysis of the observed event. These procedures rely upon sensor arrays at LIGO and Virgo to measure environmental disturbances that could potentially couple into the interferometers [67]. For all three interferometers, these procedures identified no evidence of excess power from terrestrial sources that could impact detection or analysis of GW190412. Data from Virgo contains instrumental artifacts from scattered light [68] that impact data below 20 Hz within 4 seconds of the coalescence time. As analyses of GW190412 source properties only use data from above 20 Hz, no mitigation of these artifacts was required.

III. SOURCE PROPERTIES
A. On Radiative Multipoles and Source Properties GW radiation is observed as a combination of two polarizations, h + and h × weighted with the detector re-sponse functions [69]. For GW theory, it is efficient to work with the complex valued quantity, h = h + − ih × . From the perspective of the observer, h can be expanded into multipole moments using spherical polar coordinates defined in a source centered frame [70]. Each multipole moment encodes information about the gravitationally radiating source. Interfacing GW theory with data analysis allows the connections between radiative multipole moments and source properties to be decoded.
Starting with the pioneering work of Einstein, and later refined by Newman, Penrose, Thorne and many others, GWs have been known to be at least quadrupolar [70][71][72]. The −2 spin-weighted spherical harmonics −2 Y m (θ, φ) have been found to be the simplest appropriate harmonic basis. They are the orthonormal angular eigenfunctions of gravitationally perturbed spherically symmetric spacetimes and we refer to their beyond-quadrupolar multipole moments in this basis as higher multipoles [72][73][74][75][76][77].
In this basis the multipolar decomposition is where (θ, φ) are respectively the polar and azimuthal angles defining the direction of propagation from the source to the observer, and D L the luminosity distance from the observer. The radiative multipoles, h m , depend on source properties (condensed in λ λ λ) such as the BH masses, m 1 and m 2 , and their spins S 1 and S 2 . When at least one of the BH spins is misaligned with respect to the binary's angular momentum, the orbital plane and the spins precess around the direction of the total angular momentum. We refer to these systems as precessing. For non-precessing systems, reflection symmetry about a fixed orbital plane results in a complex conjugate symmetry between moments of azimuthal index m and −m. Therefore, when we refer to a specific ( , m) multipole, h m , we always mean ( , ±m). For precessing systems and their non-fixed orbital plane, one may define a co-precessing frame such that θ is relative to the direction of maximal radiation emission. In this co-precessing frame h m approximately take on features of their non-precessing counterparts [78][79][80][81][82].
Long before merger, the instantaneous frequencies, f m , of each h m are linked to the orbital frequency of the binary f orb by f m mf orb [74]. In the moments shortly before the two BHs merge this simple scaling is broken, but its imprint persists through the final stages of coalescence [22]. Therefore, higher multipoles imply an approximately harmonic progression of frequencies within GW signals from quasi-circularly coalescing BHs.
The source geometry depends on mass ratio and is most prominently manifested in the relative contribution of multipoles with odd or even m: for an exactly equal mass binary with non-spinning components, only multipoles with even m respect orbital symmetry and so are present in the radiation [74]. In this and nearly equal mass cases, the quadrupole, h 22 , is by far the most dom-inant, followed by other multipoles with even m. However, for sufficiently unequal mass ratios, numerical and analytical studies have shown that the = m = 3 and subsequent multipoles with = m gain increasing importance [74,75,[83][84][85][86]. In the context of source detection and characterization, analytical and numerical studies [76,[87][88][89][90][91][92][93][94][95] have shown that higher multipoles can have increasing relative importance as system asymmetry increases due to source properties such as unequal masses, spin magnitudes and, in the case of precession, spin directions.
Source orientation also contributes to higher multipole content and can impact the inference of source parameters. Systems whose orbital planes point directly towards the observer present signals which are generally dominated by h 22 . In such instances the net dependence on θ and D L enter co-linearly as −2 Y 22 (θ, φ)/D L . Thus distance and inclination can be approximately degenerate for each GW detector. However, for inclined systems, the higher multipoles introduce other dependencies on θ via their harmonics, −2 Y m (θ, φ). Consequently, higher multipoles may break the inclination-distance degeneracy, thereby tightening constraints on inferred source inclination and luminosity distance. We show in Sec. III D how signal models including higher multipoles can break this degeneracy and improve our ability to infer the source parameters of GW190412. When we refer to the inclination of GW190412 below, we define it as the angle between the system's total angular momentum and the direction of propagation from the source to the observer, and denote it by θ JN for clarity.
Importantly, GW signal models are used to map between radiative multipole content and source properties, and for the robust estimation of source properties, experimental uncertainty requires GW signal models to be used in conjunction with methods of statistical inference.

B. Method and Signal Models
We perform an inference of the properties of the binary using a coherent Bayesian analysis of 8 seconds of data around the time of the detection from the three detectors (see e.g., [96]). Results presented here were mainly produced with a highly parallelized version of the code package Bilby [97,98]. Additional analyses were performed with the package LALInference [99]. RIFT [100,101] was also used to check consistency of the intrinsic parameters and for corroborating the Bayes factors that are presented below. The power spectral density (PSD) of the noise that enters the likelihood calculation is estimated from the data using BayesWave [102,103]. The low-frequency cutoff for the likelihood integration is set to 20 Hz, and the prior distributions we use are described in Appendix B.1 of [7]. We note that after initial analyses with large prior intervals on the individual masses, more restrictive prior boundaries were introduced to accelerate further Bilby analyses. Those additional bound-aries do not affect the posterior probabilities reported below as the likelihood is insignificantly small outside the allowed prior region. However, Kullback-Leibler divergences [104] calculated between prior and posterior are sensitive to this choice. Full prior specifications can be found in the data accompanying this paper [105].
The signal models we use to sample the BBH parameter space are enhanced versions of the models that have been used in past analyses (e.g., [7]). We employ models from the effective-one-body (EOB) [106][107][108][109] family that are constructed by completing an analytical inspiral-merger-ringdown description which builds on post-Newtonian (PN) [110][111][112] and black-hole perturbation theory, with numerical-relativity information. The phenomenological family [113,114] on the other hand, is based on a frequency-domain description of hybridized EOB-inspiral and numerical-relativity merger. The latest developments used here include the effects of higher multipoles in precessing models both in the EOBNR family (SEOBNRv4PHM, [33,115,116]) and the phenomenological family (IMRPhenomPv3HM, [23,24]).
All model variants that we use in the analysis of GW190412 are detailed in Table I. In order to test for imprints of spin-precession and higher multipoles in the data, we also perform analyses using models without spin precession and/or without higher multipoles. To verify the robustness of our results against waveform systematics we also performed an analysis using the numericalrelativity surrogate NRHybSur3dq8 [27] that includes the effect of higher multipoles, but is limited to spins aligned with the angular momentum. This surrogate model is constructed from numerical-relativity waveforms extended with EOB-calibrated PN waveforms.

C. Masses
In Table II we summarize the inferred values of the source parameters of GW190412. The statistical uncertainty is quantified by the equal-tailed 90% credible intervals about the median of the marginalized onedimensional posterior distribution of each parameter. We report the results obtained with the two most complete signal models -those members of the EOBNR and Phenom family that include both the effects of precession and higher multipoles (see Table I). As a conservative estimate, and because we do not favor one model over the other, we combine the posteriors of each model with equal weight, which is equivalent to marginalizing over a discrete set of signal models with a uniform probability.
The resulting values are provided in the last column of Table II, and we refer to those values in the text unless explicitly stated otherwise.
The component masses for this system in the source frame are m 1 = 29.7 +5.0 −5.3 M and m 2 = 8.4 +1.8 −1.0 M . They are consistent with the BH mass ranges of population models inferred from the first two LIGO and Virgo observing runs [7]. However, GW190412 is particularly in- The posterior distribution for the mass ratio q and effective spin χ eff of GW190412. We show the two-dimensional marginalized distribution as well as the one-dimensional marginalized distribution of each parameter along the axes for two different signal models that each include the effects of precessing spins and higher multipoles. The indicated twodimensional area as well as the horizontal and vertical lines along the axes, respectively, indicate the 90% credible regions.
teresting for its measured mass ratio of q = 0.28 +0. 13 −0.06 . Figures 2 and 3 illustrate that the mass ratio inferred for GW190412 strongly disfavors a system with comparable masses. We exclude q > 0.5 with 99% probability. In Section VI we show that the asymmetric component mass measurement is robust when analyzed using a prior informed by the already-observed BBH population.
The posteriors shown in Fig. 2 for the two precessing, higher multipole models are largely overlapping, but differences are visible. The EOBNR PHM model provides tighter constraints than Phenom PHM, and the peak of the posterior distributions are offset along a line of high correlation in the q-χ eff plane. This mass-ratio-spin degeneracy arises because inspiral GW signals can partly compensate the effect of a more asymmetric mass ratio with a higher effective spin [125][126][127][128][129]. The effective spin [108,130,131] is the mass-weighted sum of the individual spin components S 1 and S 2 perpendicular to the orbital plane, or equivalently projected along the direction of the Newtonian orbital angular momentum, L N , GW190412 is in a region of the parameter space that has not been accessed through observations before; and  we find that the two models give slightly different, yet largely consistent results. However, this is the first time that systematic model differences are not much smaller than statistical uncertainties. We tested the origin of these differences by repeating the analysis with an extended suite of signal models, as shown in Fig. 3.
The results indicate that the mass-ratio measurement of GW190412 is robust against modeling systematics, and the different treatments of higher multipoles in the EOBNR and Phenom families may account for some of the observed differences. We also see that the NRSur HM model and the EOBNR HM model agree well with each other, while the Phenom HM model deviates slightly. This is consistent with the fact that the NRSur HM and EOBNR HM models have some features in common. In NRSur HM, the PN inspiral part of the waveform is calibrated to EOB waveforms, and in EOBNR HM the merger and ringdown part of the waveform is calibrated to a subset of the numerical-relativity simulations used in the construction of NRSur HM. Further studies will be needed to fully understand the systematics visible here and mitigate them as models improve.

D. Orientation and Spins
The contribution of higher multipoles in the gravitational waveform is important for the parameter estimation of systems with small mass ratios [132,133]. In Fig. 4 we show the marginalized two-dimensional posterior distribution for luminosity distance and inclination obtained using signal models either without higher multipoles, with higher multipoles, or with higher multipoles and spin-precession. The degeneracy between luminosity distance and inclination angle that is present in the results obtained without higher multipoles is broken when higher multipoles are included. The inclusion of precession effects helps to constrain the 90% credible region further. Results obtained with the Phenom family show the same degeneracy breaking when higher multipoles are included, but the 90% credible region obtained with Phenom PHM has some remaining small support for θ JN > π/2.
We constrain the spin parameter χ eff of GW190412's source to be 0.25 +0.08 −0.11 . After GW151226 and GW170729 [2,7,34], this is the third BH binary we have identified whose GW signal shows imprints of at least one nonzero spin component, although recently another observation of a potentially spinning BH binary was reported [11]. However, inferred spins are more sensitive than other parameters (e.g., component masses) to the choice of the prior. A re-analysis of GW events with a population-informed spin prior recently suggested that previous binary component spins measurements may a Symbols: m i : individual mass; M = m 1 + m 2 ; M = (m 1 m 2 ) 3/5 M −1/5 ; superscript "det" refers to the detector-frame (redshifted) mass, while without subscript, masses are source-frame masses, assuming a standard cosmology [120] detailed in Appendix B of [7]; q = m 2 /m 1 ; M f , χ f : mass and dimensionless spin magnitude of the remnant BH, obtained through numerical-relativity fits [121][122][123][124]; χ eff , χp: effective and precessing spin parameter; χ 1 : dimensionless spin magnitude of more massive BH; D L : luminosity distance; z: redshift;θ JN : inclination angle (folded to [0, π/2]); ρ X matched-filter SNRs for the Hanford, Livingston and Virgo detectors, indicated by subscript. ρ HLV : network SNR.
have been overestimated because of the use of an uninformative prior [134]. Collecting more observations will enable us to make more confident statements on BH spins in the future. The parameter χ eff only contains information about the spin components perpendicular to the orbital plane. The in-plane spin components cause the orbital plane to precess [135], but this effect is difficult to observe, especially when the inclination angle is near 0 or π. Using models with higher multipoles, however, we constrain the inclination of GW190412 exceptionally well and put stronger constraints on the effect of precession than in previous binaries [7]. The strength of precession is parameterized by an effective precession parameter, FIG. 4. The posterior distribution for the luminosity distance, DL, and inclination, θJN (angle between the line-ofsight and total angular momentum), of GW190412. We illustrate the 90% credible regions as in Fig. 2. By comparing models that include either the dominant multipole (and no precession), higher multipoles and no precession, or higher multipoles and precession, we can see the great impact higher multipoles have on constraining the inclination and distance. All models shown here are part of the EOBNR family. 0 ≤ χ p < 1, defined by [136] . Large values of χ p correspond to strong precession. Fig. 5 shows that the marginalized one-dimensional posterior of χ p is different from its global prior distribution. The Kullback-Leibler divergence [104], D KL , for the information gained from the global prior to the posterior is 0.96 +0.03 −0.03 bits and 0.52 +0.02 −0.02 bits for the EOBNR PHM and Phenom PHM model, respectively. Those values are larger than what we found for any observation during the first two observing runs (see Table V in Appendix B of [7]). Since the prior we use introduces non-negligible correlations between mass ratio, χ eff and χ p , we check if the observed posterior is mainly derived from constraints on χ eff and q. We find that this is not the case, as a prior restricted to the 90% credible bounds of q and χ eff (also included in Fig. 5) is still significantly different from the posterior, with D KL = 0.97 +0.03 −0.03 bits (0.54 +0.02 −0.02 bits) for the EOBNR PHM (Phenom PHM) model. We constrain χ p ∈ [0.15, 0.49] at 90% probability, indicating that the signal does not contain strong imprints of precession, but FIG. 5. The posterior density of the precessing spin parameter, χp, obtained with the two models that include both the effects of precession and higher multipoles. In addition, we show the prior probability of χp for the global prior parameter space, and restricted to the 90% credible intervals of χ eff and q as given in the "Combined" column of in Table II. very small values of χ p 0.1 are also disfavored. The results obtained with the EOBNR PHM model are more constraining than the Phenom PHM results. We return to the question if GW190412 contains significant imprints of precession below, and in the context of Bayes factors in Sec. IV A.
The asymmetric masses of GW190412 means that the spin of the more massive BH dominates contributions to χ eff and χ p . Therefore, we obtain that the spin magnitude of the more massive BH is χ 1 = 0.43 +0. 16 −0.26 , which is the strongest constraint from GWs on the individual spin magnitude of a BH in a binary so far [7]. The spin magnitude of the less massive BH remains largely unconstrained.
To further explore the presence of precession in the signal, we perform the following analysis. Gravitational waveforms from precessing binaries can be decomposed into an expansion in terms of the opening angle between the total and orbital angular momenta [137,138]. Considering only = 2 modes, this expansion contains five terms, each not showing the characteristic phase and amplitude modulations of a precessing signal. When the spin component that lies in the binary's orbital plane is relatively small, the angle between the total angular momentum and the orbital angular momentum is small as well [139], and higher-order contributions in this expansion may be neglected. As a result, a precessing waveform can be modeled by the sum of the leading two contributions, where the amplitude and phase modulations of a precessing signal arise from the superposition of these terms.
In order to identify precession, we therefore require being able to measure both of these terms. We quantify the measurability of precession ρ p by how much power there is in the sub-dominant contribution. The distribution of ρ p is shown in Fig. 6. In the absence of any precession in the signal, we expect ρ 2 p to follow a χ 2 -distribution with two degrees of freedom. Using the inferred posterior distributions, our analysis shows that ρ p = 2.86 +3.43 −1.56 . We may interpret this as moderate support for precession as the median exceeds the 90% confidence interval expected from noise, but a non-negligible fraction of the ρ p posterior lies below. This calculation assumes a signal dominated by the = 2 multipole. However, we have verified that the contribution of higher multipoles to the measurement of spin-precession is subdominant by a factor of ∼5.

IV. HIGHER MULTIPOLES
Signal models that include higher multipoles are needed to infer the strongest constraints on GW190412's source properties. This is because if the data contain significant imprints of higher multipoles, the appropriate models can fit the data better than dominant-mode models, leading to a higher statistical likelihood. Conversely, if the data would not contain imprints of higher multipoles, using more complex models allows us to disfavor configurations in which clear imprints of higher multipoles are predicted [22,132,133].
In this section, we analyze how strong the imprints of higher multipoles are in GW190412 and ask if their contributions in the data are significantly stronger than random noise fluctuations. We address this question using four different approaches, each coming with its unique set of strengths and caveats.

A. Bayes Factors and Matched-Filter SNR
We may first ask if higher-multipole models actually fit the data better than dominant-multipole models. This can be quantified by the matched-filter network SNR, ρ HLV , which is based on the sum of the squared inner products between the instruments' data and the signal model. Thus the SNR quantifies the extent to which a single signal model recovers coherent power between detectors. For the EOBNR family, we find that ρ HLV increases from 18 A more complete answer to the question of which model describes the data best can be given in the Bayesian framework. The ratio of marginalized likelihoods under two competing hypotheses is called the Bayes factor, B [140]. Bayes factors may be used to quantify support for one hypothesis over another. The Bayes factor does not take into account our prior belief in the hypotheses being tested. Within GR, every compact binary coalescence signal includes higher multipoles and the prior odds in favor of their presence in the signal are infinite. We therefore focus on the Bayes factors and do not discuss the odds ratio (which is the Bayes factor multiplied by the prior odds). Table III presents log 10 B for various combinations of two models within the same waveform family. To estimate systematic uncertainties, we test the same hypotheses using multiple model families and multiple codes to calculate log 10 B. Bilby [97,98] and LALInference [99] use variants of the nested sampling algorithm [141][142][143][144]. RIFT [100,101] is based on interpolating the marginalized likelihood over a grid covering only the intrinsic source parameters. Table entries marked "-" have not been calculated because of computational constraints (LALInference analysis of precessing EOBNR models) or because the NRSurrogate model does not allow precessing spins.
We consistently find log 10 B ≥ 3 in favor of higher multipoles. This indicates strong evidence that the observed signal contains measurable imprints of higher multipoles, assuming either precessing or non-precessing (aligned spin) systems. Systematic differences between codes and waveform models dominate the uncertainty in the numbers we report. We find larger differences between codes when assuming precessing spins, because this is a more complex analysis that requires exploring more degrees of freedom in the parameter space than in the nonprecessing case. However, log 10 B is large enough across all models and codes that a statement about higher mul-tipoles can convincingly be made despite uncertainties of up to the order unity.
It would be desirable to also compare the hypotheses that the signal contains imprints of precession with assuming no precession. However, using the same codes and models that were used in Table III, the Bayes factors we found ranged from no decisive support for either hypothesis to indicating marginal support of precession. All values for log 10 B comparing precession vs. non-precessing models were smaller or comparable to the systematic uncertainties of order unity. More extensive studies will be needed to understand the origin of these systematics better.

B. Optimal SNR
A complementary way to quantify the strength of higher multipoles is to use parameter-estimation results from a signal model including higher-order multipoles [145,146]. Each multipole is decomposed into parts parallel and perpendicular to the dominant multipole by calculating the noise-weighted inner product [125,147] (often referred to as overlap) between it and the dominant multipole. Among the strongest multipoles that are included in our models, the ( , m) = (3, 3), (4,4) and (4, 3) multipoles of GW190412 are close to orthogonal to the dominant (2, 2) multipole within the band of the detector. In contrast, the (3, 2) and (2, 1) multipoles have non-negligible parallel components. To quantify the strength of the higher multipoles we remove any parallel components from the multipoles and calculate the orthogonal optimal SNR using IMRPhenomHM [22]. We find ( , m) = (3, 3) to be the strongest subdominant multipole.
The templates that include higher multipoles do not allow the amplitude and phase of the (3, 3) multipole to be free parameters; they are determined by the properties of the system. An analysis of this event using only the dominant (2, 2) multipole recovers posteriors that are consistent with a broad range of inclinations, coalescence phases, and mass ratios, while the same analysis using higher multipoles results in significantly more restricted posteriors (see Fig. 4). This suggests that by changing those parameters, our models can effectively treat the amplitude and phase of the higher multipoles as tunable parameters that make their contributions more or less pronounced. If the data only contained the dominant quadrupole mode and Gaussian noise, the squared orthogonal SNR in the subdominant multipole will be χ 2distributed with two degrees of freedom [137,138]. This was verified by analysing an injection with parameters close to GW190412.
This noise-only distribution is shown in Fig. 6, along with the orthogonal optimal SNR in the ( , m) = (3, 3) mode. The peak of the SNR distribution is at the Gaussian equivalent 3-σ level for the noise-only distribution, making this the most significant evidence for something TABLE III. log 10 B computed between two hypotheses that assume either a signal model including higher multipoles ( ≤ 5) or a dominant-multipole model ( = 2). log 10 B > 0 indicates support for higher multipoles. Each entry is based on a comparison between either precessing (first row) or nonprecessing, aligned-spin (second row) models of the same model family. See Table I for full details of the models. For each family, we also indicate the code used for calculating log 10 B. other than the dominant multipole to date [148].

C. Time-frequency Tracks
An independent analysis was performed to detect the presence of higher-order multipoles in the inspiral part of the signal, using the time-frequency spectrum of the data. Full details of the approach are described in [149], but we summarize the main idea and results for GW190412 below.
The instantaneous frequency f m (t) of the GW signal from an inspiraling compact binary is related to that of the dominant (2, 2) mode by f m (t) (m/2) f 22 (t). We compute f 22 (t) from the dominant multipoles of the EOBNR HM and Phenom HM models, using the maximum likelihood source parameters from the standard analysis presented in Sec. III. Inspired by the above scaling relation, we then look along the generalized tracks, f α (t) αf 22 (t), in a time-frequency representation that is the absolute square of the continuous wavelet transformation (CWT) of the whitened on-source data,X(t, f ). We have used Morlet wavelets to perform the CWT, where the central frequency of the wavelet was chosen so as to maximize the sum of pixel values along the f 22 (t) curve. This wavelet transformation is shown in the top panel of Fig. 7.
In order to quantify the energy along each track, we define Y (α) to be the energy |X(t, f )| 2 in all pixels containing the track f α (t), where we discretize the data with a pixel size of ∆t = 1/4096 s along the time-axis and ∆f = 1/5 Hz along the frequency axis. By doing so, we can decouple the energy in individual multipoles of the signal. Once f 22 (t) is defined, this is a computationally efficient way to analyze which multipoles have sufficient energy to be detectable in the data; no further modeling input is needed, although we do not require phase coherence along each track.
The resulting Y (α) for GW190412 is shown in the bottom panel of Fig. 7. It has a global peak at α = 1, corresponding to the dominant (2, 2) multipole, and a prominent local peak at α = 1.5, corresponding to the m = 3 multipoles. We also calculate Y (α) from different segments surrounding, but not including GW190412, to capture the detector noise characteristics; in this case we call the quantity N (α). The ensemble average µ(α) = N (α) and standard deviation σ(α) of N (α) are also plotted for reference and to highlight the relative strength of the GW signal present in the on-source segment.
Instead of estimating the significance of the m = 3 multipoles from comparing Y (α) to its background at α = 1.5, we perform a more powerful statistical analysis in which we test the hypotheses that the data contain either only noise (H 0 ), or noise and a dominant-multipole maximum likelihood signal (H 1 ), or noise and a maximum likelihood signal that includes m = 2 and m = 3 multipoles (H 2 ). By maximizing the likelihood of observing Y (α) given each hypothesis over a free amplitude parameter for each multipole, we obtain likelihood ratios for H 1 and H 2 , and their difference is in turn incorporated into a detection statistic β (see Eq. (7) in [149]).
From the on-source data segment taken from the LIGO Livingston detector, we found the detection statistic β = 6.5 with a p-value of 6.4 × 10 −4 for EOBNR HM model; and β = 6.6 with a p-value < 6.4 × 10 −4 for Phenom HM model, which strongly supports the presence of m = 3 modes in the signal. The full distribution of β from offsource data segments from the LIGO Livingston detector surrounding the trigger time of GW190412 is shown in the inset of Fig. 7.

D. Signal Reconstructions
As an additional test of consistency, and an instructive visual representation of the observed GW signal, we compare the results of two signal reconstruction methods. One is derived from the parameter-estimation analysis presented in Sec. III, the second uses the model-agnostic wavelet-based burst analysis BayesWave [150] which was also used to generate PSDs. A detailed discussion of such signal comparisons for previous BBH observations can be found in [151].
For GW190412, both signal reconstruction methods agree reasonably well as illustrated in Fig. 8. To quantify the agreement for each signal model from the Phenom and EOB families, we compute the noise-weighted inner product [125,147] between 200 parameter-estimation samples and the BayesWave median waveform. The BayesWave waveform is constructed by computing the median values at every time step across samples. Similar comparison strategies have been used in [3,5,34,152]. For the models of the EOBNR family, we find that the agreement with the unmodeled BayesWave reconstruction increases slightly from overlaps of 0.84 +0.01 −0.02 for the dominant-multipole, non-precessing model to 0.86 +0.01 −0.02 when higher multipoles are included, to 0.88 +0.01 −0.02 for the most complete EOBNR PHM model (median overlap and 90% errors). Increasing overlaps are consistent with the findings of the other methods presented in this section that indicate that the extra physical effects included in the higher-multipole precessing model match the data better. The overlaps we find are also consistent with expectations from [151]. The Phenom family may suggest a similar trend. The overlap for non-precessing dominant multipole model is 0.84 +0.02 −0.02 , and it increases to 0.85 +0.01 −0.02 for the non-prcessing higher multipoles model, to 0.86 +0.02 −0.03 for the precessing higher multipoles model Phenom PHM.

V. TESTS OF GENERAL RELATIVITY
As the first detected BBH signal with a mass ratio significantly different from unity, GW190412 provides the opportunity to test GR in a previously unexplored regime. Due to the mass asymmetry, this signal contains information about the odd parity multipole moments. Hence the tests of GR reported here are sensitive to potential deviations of the multipolar structure away from GR [153]. A violation from GR may arise from how the signal is generated by the source; additionally, the form of the signal described by GR may be tested by checking the consistency of independently obtained estimates of parameters between the inspiral and merger-ringdown parts of the full BBH waveform. The following analyses are done by using the LALInference library [99] to generate posterior probability distributions on these parameters by using the nested sampling algorithm.

A. Constraints on gravitational wave generation
We check the consistency of this source with general relativistic source dynamics by allowing for parameterized deformations in each phasing coefficient in the binary's waveform. They were first performed on inspiralonly waveforms in [154,155] and an extension to highermodes was studied in [156]. The current version of the test using the phenomelogical waveform models rely on extensions of the methods laid out in [157,158]. Such tests have been performed on all GW detections made in O1 and O2 [159,160] and have been updated recently with the best constraints by combining all significant BBH detections made during O1 and O2 in [161]. We perform this analysis with the precessing, dominantmultipole phenomenological model Phenom P, and, for consistency with previous tests, with the aligned-spin dominant-mode EOBNR model (see Table I). . We show detector data, whitened by an inverse amplitude-spectral-density filter computed using BayesLine [102], together with the unmodeled BayesWave reconstruction that uses a wavelet bases, and the reconstruction based on the precessing, higher multipole models from the EOBNR and Phenom families. The bands indicate the 90% credible intervals at each time. We caution that some apparent amplitude fluctuations in this figure are an artifact of the whitening procedure.
The inspiral regime of both waveforms is modeled using the PN approximation. The fractional deviation parameters δφ n are added to their respective PN coefficients ϕ n at n/2-PN order. While the deviation coefficients are added differently in the two models, the differences are taken care of by effectively re-parameterizing the coefficients added to the EOB-based model for consistency in comparing bounds from the Phenom-based model. The full set of parameters being tested are {δφ 0 , δφ 1 , δφ 2 , δφ 3 , δφ 4 , δφ 5l , δφ 6 , δφ 6l , δφ 7 }. Here, δφ 5l and δφ 6l refer to the fractional deviations added to the log-dependent terms at 2.5PN and 3PN respectively. Moreover, δφ 1 is an absolute deviation as there is no 0.5PN term within GR. These parameters are tested by varying only one δφ n at a time, and introducing these to the parameter set of the full signal model. This increased parameter space dimensionality makes it especially challenging to use the already computationally expensive Phenom PHM waveform, and we restrict our analysis to using the precessing Phenom P approximant. To check for possible systematics introduced from analyzing this signal with the dominant-mode Phenom P model, a signal similar to GW190412 was created with the Phenom PHM waveform model and injected into data generated using the BayesWave PSDs for the event. The recovery of this signal, using Phenom P shows that the posteriors are completely consistent with the injected values, suggesting that for a GW190412-like signal with same SNR, the inclusion of the higher multipoles does not significantly bias results when those multipoles are not included.
The posterior distributions on the fractional deviation parameters are always found to be consistent with the GR prediction of δφ n = 0. Additionally, the EOB model can test deviations in the −1PN dipole term, while the phenomenological signal model can be used to test the intermediate and merger-ringdown parameters of the signal, which are consistent with their GR value. How- ever, owing to the longer inspiral of this signal, the bounds obtained from this event in the inspiral regime are among the most constraining bounds obtained from the analyses on individual BBH detections as reported in [161]. We show the 90% upper bounds computed in Fig. 9. The only BBH signals that give more robust or comparable bounds above 0 PN order are those from GW170608 (the lowest mass BBH to have been published) and GW151226.

B. Inspiral-Merger-Ringdown Consistency
We check the consistency of signal parameters between the low and high frequency parts of the signal [162,163].
Estimates of the final mass M f and final spin χ f of the remnant BH are found from the two parts of the frequency domain signal and their fractional differences are checked for consistency. For this source, the transition from the lower-frequency to higher-frequency part (the Kerr innermost stable circular orbit) of the signal is estimated from the median intrinsic source parameters and the resulting prediction for M f and χ f to be at f = 211 Hz [164]. We used the signal model Phenom PHM, sampling on the BBH parameter set to obtain posterior probability distributions on all parameters. The component masses and spins are estimated directly from the lower-frequency part of the signal, and, using numerical relativity fit formulae [121][122][123][124], the posteriors on M f and χ f are inferred. From the higher-frequency part of the signal, estimates on component masses and spins are obtained again using the same waveform model, and the posteriors on M f and χ f are inferred using the same fit formulae as above. From those two distributions, a posterior distribution of the fractional differences, denoted by ∆M f /M f and ∆χ f /χ f respectively, is then computed. Here,M f andχ f denote the mean values of the distributions of M f and χ f respectively. While we expect mass and spin differences of exactly (0,0) given a pure signal obeying GR, the presence of detector noise will generally yield a posterior with some non-zero spread around (0,0). Figure 10 shows the results of the posterior distributions on these fractional quantities. The 90% credible regions of the quantities ∆M f /M f and ∆χ f /χ f enclose (0,0) as can be seen from both the one-dimensional posteriors and the two-dimensional contours. GW190412 results are consistent with past observations of BBHs [161].

VI. IMPLICATIONS FOR BBH POPULATION PROPERTIES
BBHs detected by the LIGO-Virgo network can be used to constrain the uncertain physical processes inherent to compact binary formation channels. As the first observed BBH with definitively asymmetric masses, the inclusion of GW190412 in the current catalog of BBHs has a significant impact on inferred population properties. Here, we examine (i) how the addition of GW190412 to the catalog of BBHs from the first and second observing runs affects population statements; (ii) the robustness of the component mass measurements of GW190412 when evaluated as part of the previously observed population; and (iii) whether this system's mass ratio is a significant outlier with respect to that population.
Using the ten significant BBH events in the catalog of GWs from the first and second observing runs (GWTC-1, [7]), we constrained the parameters of phenomenological models that represent the underlying BBH population [15]. In certain models, the mass-ratio distribution is parameterized with a power law, p(m 2 |m 1 ) ∝ q βq [165][166][167], where β q is the spectral index of the mass ratio dis- tribution. Since all 10 events from GWTC-1 are consistent with symmetric masses, the posteriors for β q showed a preference for positive values [15], providing initial evidence supporting equal-mass pairings over randomlydrawn mass pairings [168]. However, the steepness of β q was unconstrained, which limited our ability to determine how prevalent equal-mass pairings are relative to their asymmetric counterparts. The inclusion of GW190412 in the population provides a much stronger constraint on the mass ratio spectral index, as shown in Figure 11. Applying the population of significant events from the first and second observing runs as well as GW190412 to the simplest model that invokes a power-law spectrum to the mass ratio distribution (Model B from [15]), we find β q < 2.7 (5.8) at the 90% (99%) credible level. This indicates that though equal-mass pairings may still be preferred, there is significant support for asymmetric mass pairings; the posterior population distribution [15,169] for this model indicates that 10% of merging BBHs should have a mass ratio of q 0.40. In fact, the support for β q ≤ 0 in the recovered distribution indicates that the true mass ratio distribution may be flat or even favor unequal-mass pairings. This is not in tension with the mass ratios of the already-observed population; though all mass ratio posteriors from GWTC-1 are consistent with q = 1, they also have significant support for lower values. These constraints on β q are preliminary and final results from O3 will only be obtained after analyzing the population that includes all BBH events from this observing run.
FIG. 11. Posterior on mass ratio spectral index βq with (solid lines) and without (dashed lines) the inclusion of GW190412. We show inference using both the EOBNR and Phenom families; for the 10 events from GWTC-1 we use the publicallyavailable samples for both of these waveform families presented in [7], and for GW190412 we use the EOBNR PHM and Phenom PHM posterior samples presented in this paper.
We also check whether the asymmetric mass ratio for GW190412 is robust when the component mass posterior distributions are reweighted using a population-informed prior based on Model B from [15]. Since the majority of observed systems are consistent with equal mass, the mass ratio posterior for GW190412 pushes closer towards equal mass when using a population-informed prior rather than the standard uninformative priors used to generate posterior samples. However, the mass ratio of GW190412 is still constrained to be q < 0.43 (0.59) at the 90% (99%) credible level, compared to q < 0.38 (0.49) using the standard priors from parameter estimation.
Using methods from [170], we test the consistency of GW190412 with the population of BBHs inferred from the first and second observing runs. We first construct a population model (Model B) derived only from the events in GWTC-1, following the prescription in [15], and draw 11 observations from this model (representing the 10 significant BBHs from the first two observing runs as well as GW190412). Examining the lowest mass ratio drawn from each set of 10 7 such realizations, we find the population-weighted mass ratio posterior samples of GW190412 lie at the 1.7 +10.3 −1.3 percentile of the cumulative distribution of lowest mass ratios. This indicates that given the BBH population properties inferred from the first two observing runs, drawing a system with a mass ratio analogous to GW190412 would be relatively rare. The apparent extremity of GW190412 is likely driven by the limited observational sample of BBHs and the lack of constraining power on the mass ratio spectrum prior to the observation of GW190412. Constructing a population model that includes the observation of GW190412 in the fit, we find GW190412 to lie at the 25 +47 −17 percentile of the cumulative distribution of lowest mass ratios, in-dicating that GW190412 is a reasonable draw from the updated population.

VII. ASTROPHYSICAL FORMATION CHANNELS FOR GW190412
Multiple astrophysical channels are predicted to produce the merging BBHs identified by the LIGO-Virgo network. The majority of these channels have mass ratio distributions that peak near unity, but also often predict a broad tail in the distribution that extends to more asymmetric masses. Though a wide array of formation channels exist, each with distinct predictions for merger rates and distributions of intrinsic BBH parameters, most channels can be broadly categorized as the outcome of either isolated binary stellar evolution or dynamical assembly (for reviews, see [171,172], respectively).
In the canonical isolated binary evolution scenario, by which a compact binary progenitor achieves a tight orbital configuration via a common envelope phase [171,[173][174][175][176][177][178], BBH mergers with mass ratios of q 0.5 are typically found to be less common than their nearequal-mass counterparts by an order of magnitude or more [179][180][181][182][183][184], though certain population synthesis modeling finds BBH mergers with asymmetric masses to be more prevalent [185,186]. However, even if the formation probability of asymmetric mass ratio systems is an order of magnitude lower than the formation probability of equal mass systems, the observation of one clearly asymmetric mass system given the current observational catalog is unsurprising. In contrast, the progenitor of GW190412 is unlikely to have formed through chemically homogeneous evolution, as this scenario typically cannot form binaries with mass ratios below q < 0.5 [187][188][189]. The asymmetric component masses of GW190412 may also suggest formation in an environment with lower metallicity, as lower metallicities are predicted to produce a higher rate of merging BBHs having mass ratios consistent with GW190412 [178,179,183,184], though this prediction is not ubiquitous across population synthesis models [190].
Dense stellar environments, such as globular clusters [191][192][193][194], nuclear clusters [195], and young open star clusters [196][197][198][199][200], are also predicted to facilitate stellar-mass BBH mergers. Numerical modeling of dense clusters suggests that significantly asymmetric component masses are strongly disfavored for mergers involving two first-generation BHs that have not undergone prior BBH mergers (e.g., [201]). However, asymmetric component masses may be attained by a first-generation BH merging with a higher-generation BH that has already undergone a prior merger or mergers [166,202,203]. For formation environments such as globular clusters, this would require low natal spins for the initial population of BHs so that an appreciable number of merger products can be retained in the shallow gravitational potential of the cluster [203,204]. Alternatively, BBH mergers with asymmetric masses may be the result of massive-star collisions in young star clusters, as these environments have been shown to amplify unequal-mass BBH mergers relative to their isolated counterparts [199].
Other formation scenarios may also be efficient at generating BBH mergers with significantly asymmetric component masses. BBHs in triple or quadruple systems can undergo Lidov-Kozai oscillations [205,206] that may expedite the GW inspiral of the binary. Such systems can either exist in the galactic field with a stellar-mass outer perturber [207][208][209][210][211][212][213] or in galactic nuclei with a supermassive BH as the tertiary component [214][215][216][217]. Though most modeling of hierarchical stellar systems do not include robust predictions for mass ratio distributions of merging BBHs, certain models find galactic field triples with asymmetric masses for the inner BBH to have a merger fraction that is about twice as large as their equalmass counterparts [207]. In the context of hierarchical triples in galactic nuclei, recent modeling predicts a significant tail in the mass ratio distribution of merging BBHs that extends out to mass ratios of ∼ 8:1 [217]. BBH mergers with significantly asymmetric component masses are also predicted for systems formed in the disks of active galactic nuclei [218][219][220][221][222][223][224]. The deep gravitational potential near the vicinity of the supermassive BH may allow for stellar-mass BHs to go through many successive mergers without being ejected by the relativistic recoil kick, leading to BBH mergers with highly asymmetric masses [223,224]. Though these channels may not be dominant, they could explain a fraction of the sources observed by the LIGO-Virgo network.
In summary, though the mass ratio of GW190412 is the most extreme of any BBH observed to date, it is consistent with the mass ratios predicted from a number of proposed BBH formation channels. Many astrophysical channels predict that near-equal-mass BBH mergers are more common than mergers with significantly asymmetric component masses. However, as the observational sample of BBHs grows, it is not unexpected that we would observe a system such as GW190412 that occupies a less probable region of intrinsic parameter space. Future detections of BBHs will enable tighter constraints on the rate of GW190412-like systems.

VIII. CONCLUSIONS
Every observing run in the advanced GW detector era has delivered new science. After the first observations of BBHs in the first observing run, and the continued observation of BBHs as well as the multimessenger observation of a binary neutron star in the second observing run [6,225], O3 has been digging deeper into the populations of compact binary mergers. The observation of a likely second neutron-star binary in O3 has already been published [8], and here we presented another GW observation with previously unobserved features.
GW190412 was a highly significant event, with a com-bined SNR of 19 across all three GW detectors. It is the first binary observed that consists of two BHs of significantly asymmetric component masses. With 99% probability, the primary BH has more than twice the mass of its lighter companion. The measurement of asymmetric masses is also robust even when the properties of GW190412 are inferred using a population-based prior. This observation indicates that the astrophysical BBH population includes unequal-mass systems. GW190412 is also a rich source from a more fundamental point of view. GR dictates that gravitational radiation from compact binaries is dominated by a quadrupolar structure, but it also contains weaker contributions from subdominant multipoles. Here we provided conclusive evidence that at least the second most important multipole -the ( , |m|) = (3, 3) multipole -makes a significant, measurable contribution to the observed data. As a result, the orientation of the binary is more accurately determined and tighter bounds are obtained on relevant intrinsic source parameters such as the mass ratio and spin of the system.
The asymmetric mass ratio of GW190412 allows the primary spin to have a more measurable effect on the signal. We find the primary spin magnitude to be 0.43 +0. 16 −0.26 , which is the strongest constraint on the individual spin magnitude of a BH using GWs so far. Though we only find marginal statistical hints of precession in the data, the results presented here illustrate that we confidently disfavor strong precession (as would be characterized by a large in-plane spin parameter).
GW190412 is a BBH that occupies a previously unobserved region of parameter space. As we continue to increase the sensitivity of our detectors and the time spent observing, we will gain a more complete picture of the BBH population. Future observations of similar types of binaries, or even more extreme mass ratios, will sharpen our understanding of their abundance and might help constrain formation mechanisms for such systems. GW190412 also shows that numerical and analytical advances in modeling coalescing binaries in previously unexplored regimes remains crucial for the analysis of current and future GW data. The most recent and most complete signal models robustly identified GW190412's source properties and showed consistency with GR. Systematic differences are visible and will become more important when we observe stronger signals, pointing to the necessity for future work in this area.
LIGO and Virgo data containing GW190412, and samples from a subset of the posterior probability distributions of the source parameters [105], are available from the Gravitational Wave Open Science Center [226].
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 Foundation for Fundamental Research on Matter supported by the Netherlands Organisation 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