Observation of Gravitational Waves from Two Neutron Star-Black Hole Coalescences

We report the observation of gravitational waves from two compact binary coalescences in LIGO ’ s and Virgo ’ s third observing run with properties consistent with neutron star – black hole ( NSBH ) binaries. The two events are named GW200105_162426 and GW200115_042309, abbreviated as GW200105 and GW200115; the ﬁ rst was observed by LIGO Livingston and Virgo and the second by all three LIGO – Virgo detectors. The source of GW200105 has component masses -+ 8.9 1.51.2 and  -+ M 1.9 0.20.3 , whereas the source of GW200115 has component masses -+ 5.7 2.11.8 and  -+ M 1.5 0.30.7 ( all measurements quoted at the 90% credible level ) . The probability that the secondary ’ s mass is below the maximal mass of a neutron star is 89% – 96% and 87% – 98%, respectively, for GW200105 and GW200115, with the ranges arising from different astrophysical assumptions. The source luminosity distances are

8.9 1.5 1.2 and  -+ M 1.9 0.2 0.3 , whereas the source of GW200115 has component masses (all measurements quoted at the 90% credible level).The probability that the secondary's mass is below the maximal mass of a neutron star is 89%-96% and 87%-98%, respectively, for GW200105 and GW200115, with the ranges arising from different astrophysical assumptions.The source luminosity distances are , respectively.The magnitude of the primary spin of GW200105 is less than 0.23 at the 90% credible level, and its orientation is unconstrained.For GW200115, the primary spin has a negative spin projection onto the orbital angular momentum at 88% probability.We are unable to constrain the spin or tidal deformation of the secondary component for either event.We infer an NSBH merger rate density of - + -- 45 Gpc yr 33 75 3 1 when assuming that GW200105 and GW200115 are representative of the NSBH population or -+ --

Introduction
In 2020 January, the LIGO-Virgo detector network observed gravitational-wave (GW) signals from two compact binary inspirals that are consistent with neutron star-black hole (NSBH) binaries.These represent the first confident observations to date of NSBH binaries via any observational means.The two events, carrying the full designations GW200105_162426 and GW200115_042309 and abbreviated henceforth as GW200105 and GW200115, were detected on 2020 January 5 at 16:24:26 UTC and 2020 January 15 at 04:23:10 UTC, respectively.The coincident detection of GW200115 by the three detectors LIGO Hanford, LIGO Livingston, and Virgo gives it high confidence of being an astrophysical GW event.During the other event, GW200105, the LIGO Hanford detector was not operational, and, owing to the small signal-to-noise ratio (S/N) in Virgo, it was effectively a single-detector event in LIGO Livingston.While quantification of the confidence of single-detector events is subject to significant uncertainty, GW200105 stands clearly apart in the LIGO Livingston data from any other candidate with NSBH-like parameters during the ∼11 months of the third observing run (O3).
The component masses inferred from these binary inspirals provide information on the nature of their components.The primaries have masses of = -+ m 8.9 1 1.5 1.2 and  -+ M 5.7 2.1 1.8 for GW200105 and GW200115, respectively, with uncertainties quoted at the 90% credible level.These primary masses are well above the maximum mass of a neutron star (NS; Rhoades & Ruffini 1974;Abbott et al. 2018a;Cromartie et al. 2019;Shibata et al. 2019;Farr & Chatziioannou 2020;Fonseca et al. 2021;Nathanail et al. 2021) and within the mass range of black holes (BHs) observed electromagnetically (Özel et al. 2010; Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence.Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.Farr et al. 2011;Kreidberg et al. 2012;Miller & Miller 2014) and via GWs (Abbott et al. 2016a(Abbott et al. , 2019c(Abbott et al. , 2021b)) , respectively, for GW200105 and GW200115, within the mass range of known NSs (Antoniadis et al. 2016;Abbott et al. 2017aAbbott et al. , 2020a;;Alsing et al. 2018).
Detections of NSBHs have so far remained elusive in both electromagnetic (EM) and GW surveys.In the past four decades, surveys have identified 19 binary neutron star (BNS) systems in the Galaxy (Farrow et al. 2019;Agazie et al. 2021); however, the discovery of a pulsar in an NSBH binary remains a key objective for current and future radio observations (Liu et al. 2014;Weltman et al. 2020).
Similarly, GW observations of LIGO and Virgo through the first part of the O3 run (O3a) have led to the identification of 48 binary black hole (BBH) candidates (Abbott et al. 2019c(Abbott et al. , 2021b) ) and two BNS candidates (Abbott et al. 2017a(Abbott et al. , 2020d)).Independent analyses of the public detector data identified additional GW candidates (Magee et al. 2019;Nitz et al. 2019bNitz et al. , 2020aNitz et al. , 2021;;Venumadhav et al. 2019Venumadhav et al. , 2020;;Zackay et al. 2019aZackay et al. , 2019b)).The absence of NSBH candidates in LIGO's and Virgo's first and second observing runs (O1 and O2, respectively) led to the upper limit on the local merger rate density of NSBH systems of   610 NSBH Gpc −3 yr −1 at the 90% credible level (Abbott et al. 2019c).
During O3a, two events were notable as possible NSBH candidates.First, GW190426_152155 (Abbott et al. 2021b) was identified as a marginal NSBH candidate with a falsealarm rate (FAR; 1.4 yr −1 ) so high that it could also plausibly be a detector noise artifact.Second, GW190814 (Abbott et al. 2020c) (Essick & Landry 2020;Fattoyev et al. 2020;Godzieba et al. 2021;Tews et al. 2021), such as those found in known binaries that will coalesce within a Hubble time, the secondary could conceivably be an NS spinning near its breakup frequency (Essick & Landry 2020;Most et al. 2020;Dexheimer et al. 2021;Tews et al. 2021).
This Letter presents the status of the detectors during times around GW200105 and GW200115 (Section 2) and the results of the different searches leading to the detections (Section 3).We describe the main properties of the two events (Section 4) and discuss the nature of the secondary components (Section 5).Finally, we present the astrophysical implications, including the merger rates of this new class of GW source (Section 6), and conclude (Section 7).Data products associated with the events reported here, such as calibrated strain time series and parameter estimation posterior samples, are available through the Gravitational Wave Open Science Center (GWOSC) at gw-openscience.org.

Detectors and Data
The two events reported here were observed in the second part of the third observing run (O3b).At the time of GW200105, LIGO Livingston had been in a stable operational state for over 10 hr, with a sensitivity, quantified by the angleaveraged BNS inspiral range (Allen et al. 2012), of ∼137 Mpc.Virgo had been in its nominal state for ∼22 hr with a BNS range of ∼45 Mpc, and LIGO Hanford was not operational.At the time of GW200115, all three interferometers had been in a stable operational state for over 2 hr.The BNS ranges for LIGO Hanford, LIGO Livingston, and Virgo around the detection were ∼115, ∼133, and ∼50 Mpc, respectively.
Figure 1 shows time-frequency representations (Chatterji et al. 2004) of the two events.The LIGO Livingston data of GW200105 show a track of excess power with increasing frequency.For GW200115, no similar tracks are visible, as the S/N in each of the detectors is lower than that of GW200105 in LIGO Livingston (see Figure 3).Also, light-scattering noise (Soni et al. 2020(Soni et al. , 2021) ) is visible in Figure 1 in LIGO Livingston around 20 Hz.
The LIGO and Virgo GW detectors are calibrated using radiation pressure from auxiliary lasers at a known frequency and amplitude (Acernese et al. 2018;Sun et al. 2020).In LIGO, the calibration pipeline (Viets et al. 2018) linearly subtracts the noise from the calibration lines and the harmonics of the power mains from the data.
The calibration systematic error and associated uncertainty of the data in the [20-1024] Hz frequency region are for the LIGO detectors no larger than 8.6% in amplitude and 5.9°in phase for GW200105 and 8.0% in amplitude and 5.5°in phase for GW200115 (68% credible intervals).For Virgo, we use 5.0% in amplitude and 4.1°in phase for both events, except for frequencies 46-51 Hz, where additional calibration systematic error arises from online loops set up to damp the main power lines at 50 Hz and mechanical resonances of the suspensions close to 48 Hz.This effect was introduced during detector improvements carried out between O3a and O3b and is not accounted for in the signal reconstruction process.However, Virgo data in the frequency window 46-51 Hz are excluded from the parameter estimation analyses in Section 4.
To verify that instrumental noise artifacts do not bias the analysis of source properties of the observed events, we use data quality validation procedures as in previous events (Abbott et al. 2016e;Davis et al. 2021), employing sensor arrays at LIGO and Virgo to measure environmental disturbances that could couple into the interferometers (Nguyen et al. 2021).In Virgo, we find no evidence of excess power from terrestrial sources for both events.For GW200105, we identify lightscattering noise in LIGO Livingston below 25 Hz, 3 s before merger.As in past detections (Abbott et al. 2021b), we subtract this noise with the BayesWave algorithm (Cornish et al. 2021) and use the cleaned data for the source parameter estimation in Section 4. The top left panel of Figure 1 shows the cleaned data.A low-energy feature around 20 Hz, not overlapping with the time-frequency track of GW200105, remains after the glitch subtraction.For GW200115, we also identify light scattering in the LIGO Livingston data below 25 Hz.Due to the increased difficulty of subtracting the long-duration glitching that coincides with the time-frequency track of GW200115, glitch-subtracted data were not available at the time of analysis.Hence, we exclude LIGO Livingston data below 25 Hz in the analysis in Section 4.

Detections
Event GW200115 is a coincident event, and in Section 3.1, we describe the established procedures to identify it and determine its significance.Subsequently, we describe the procedures followed for the single-detector event GW200105.
A public Preliminary Gamma-ray burst Coordinates Network (GCN) Notice was sent out 6 minutes after the event (LIGO Scientific Collaboration & Virgo Collaboration 2020a).The lowlatency BAYESTAR (Singer & Price 2016) sky map computed from the original trigger, which was produced from MBTA ONLINE, indicated a possible discrepancy compared to those from the other three pipelines.Therefore, the representative trigger was manually switched to the one from GSTLAL (LIGO Scientific Collaboration & Virgo Collaboration 2020b).The associated sky map is shown in the bottom panel of Figure 2 with a solid black line with 90% credible area of 900 deg 2 .The GCN Circular reported a probability >99% for the lighter compact object to have a mass below 3 M e , derived from a machine-learning analysis (Chatterjee et al. 2020;Kapadia et al. 2020).The RAVEN pipeline (Urban 2016), which looks for coincident gamma-ray bursts (GRBs), did not report an associated GRB trigger from Fermi-GBM or the Neil Gehrels Swift Observatory within −1 and +5 s of the GW trigger.A total of 31 GCN Circulars (GCN archive for S200115j; 2020) reporting follow-up observations were published for this event.To date, no EM counterpart has been reported associated with GW200115.
After the low-latency identification of GW200115, the extended periods of strain data around the event are further analyzed by the detection pipelines in their offline configurations using improved calibration and refined data quality information that is not available in low latency (Abbott et al. 2018b(Abbott et al. , 2021b;;Davis et al. 2021).In this analysis, we also subtract nonstationary noise due to the nonlinear coupling of the 60 Hz power mains at the LIGO detectors using coupling functions derived from machine-learning techniques (Vajente et al. 2020).
The PYCBC pipeline uses a template bank constructed with a hybrid geometric-random algorithm (Roy et al. 2017(Roy et al. , 2019)), while GSTLAL and MBTA use a template bank generated by a stochastic placement method (Babak 2008;Privitera et al. 2014;Mukherjee et al. 2021).The GSTLAL pipeline identifies GW200115 with a network S/N of 11.6 and FAR of <1/(1 × 10 5 yr) using the data from 2019 November 1 15:00 UTC to 2020 January 22 18:11 UTC.The MBTA and PYCBC offline analyses also identified the trigger in data from 2020 January 13 10:27 UTC to 2020 January 22 18:11 UTC, yielding consistent network S/Ns as shown in Table 1, as well as estimated FARs of 1/(182 yr) and <1/(5.6 × 10 4 yr), respectively.Based on two detectors, these FARs robustly indicate the confidence of the detection.

GW200105-Single-detector Event
Event GW200105 was identified by GSTLAL running in the low-latency configuration as a possible CBC candidate in LIGO Livingston and Virgo data with a network S/N of 13.9 and merger time 2020 January 5 at 16:24:26 UTC.The candidate identified as S200105ae in GRACEDB was considered a single-detector event because the S/N in Virgo was below the threshold S/N of 4.0.
The low-latency FAR of ≈24/yr did not pass the threshold for sending a GCN alert.Without the information provided by temporal coincidence between different detectors and consistency in recovered parameters, it is more difficult to obtain a robust estimate of the event's significance.In fact, an alternative configuration of GSTLAL, running to test several improvements made in the offline configuration of early O3 (Abbott et al. 2021b), found the trigger at higher significance.Therefore, validation procedures and subsequent early-offline analysis were initiated.The three other low-latency search pipelines, MBTA ONLINE, PYCBC LIVE, and SPIIR, also generated triggers with consistent S/Ns near the event time of S200105ae (see Table 1).These three pipelines were not configured to assign FARs to single-detector triggers and therefore did not generate automatic alerts for GW200105.
On 2020 January 6 at 19:39:09 UTC, the initial GCN Circular reported a low-latency BAYESTAR skymap (solid black contours) in the top panel of Figure 2 with 90% credible region of 7700 deg 2 and a 12% chance of formation of tidally disrupted matter-relevant for a possible EM counterpart-as a result of the coalescence, as derived from a machine-learning analysis (Chatterjee et al. 2020;Kapadia et al. 2020).Using the parameters obtained from low-latency parameter estimation, the chance of formation of tidally disrupted material was later revised to be negligible (see LIGO Scientific Collaboration & Virgo Collaboration 2020c and Section 5 below).
In response to the discovery notice of GW200105, multiple observatories carried out follow-up and shared their results publicly via 21 GCN circulars (GCN archive for S200105ae 2020).To date, no EM counterpart has been reported associated with GW200105.
After the event's identification in low latency, the strain data are further analyzed by the detection pipelines in their offline configuration, following procedures for single-detector events discussed in Abbott et al. (2020a).Here GSTLAL identifies GW200105 with an S/N of 13.9 and FAR of 1/(2.8 yr) in its offline configuration using data from 2019 November 1 15:00 UTC up to 2020 January 22 18:11 UTC.For single-detector events, like GW200105, the FAR estimate involves extrapolation, as explained below.The MBTA and PYCBC offline analyses also identify the trigger in the LIGO Livingston data from 2020 January 4 17:09 UTC to 2020 January 13 10:27 UTC yielding the consistent S/Ns shown in Table 1.The S/N for GW200105 is larger than that for GW200115, even though LIGO Hanford was not operational.
We cannot rely on only the S/N of the CBC trigger to estimate its significance, i.e., to quantify how often detector noise mimics a possible coalescence signal.Because detector noise shows significant deviation in the tails from the standard Gaussianity assumptions, its properties have to be estimated empirically (Abbott et al. 2016b).For multi-detector triggers like GW200115, consistency of the observed CBC trigger among two or more detectors forms an integral part of establishing it as a confident detection with a low FAR.On the other hand, single-detector triggers can be assigned only modest FAR values due to limited detector observational time (Callister et al. 2017;Nitz et al. 2020b).Nonetheless, alternative methods can be used to estimate confidence in a single-detector event, as was demonstrated for GW190425 (Abbott et al. 2020a).
The signal GW200105 measured in LIGO Livingston is distinct from all noise events in the entire O3 observation period.Specifically, Figure 3 shows with colored shading the distribution of background noise triggers identified by GSTLAL in the region of binary systems with a chirp mass less than 4 M e .The colored region indicates the density of background noise triggers quantified through the S/N-ξ 2 noise probability density function for the three detectors, where the autocorrelation ξ 2 provides a consistency test similar to χ 2 (Messick et al. 2017).The red star represents GW200105, which lies in a region of this plane without noise triggers, indicating that this trigger is unique within the entire O3.
For comparison, Figure 3 also shows GW200115 and the marginal candidate GW190426_152155 separately for LIGO Livingston and LIGO Hanford.In general, triggers such as GW200115 and GW190426_152155 in Figure 3 would not be separable from the noise by a single detector alone, and it is the coincidence between multiple detectors that raises their significance.On the other hand, the coincidence was not available for GW200105, yet it clearly stands out from the background in Figure 3, as strong GW signals typically do.
The uniqueness of GW200105 within the data shown leads to a upper bound on the FAR assignment comparable to the inverse observing time.A stronger bound, like the FAR quoted above, must rely on assumptions made on the properties of the noise triggers if one had collected the observational strain data for a longer period, a process called extrapolation.As a consequence of the assumptions made, the extrapolated FAR estimates have large uncertainties.Since single-detector FAR estimates rely on the uniqueness within the data set, the values may change and errors may reduce as more data are accumulated.At the time of the analysis, neither MBTA nor PYCBC had the capability to assign a significance to singledetector triggers, although GW200105 also stands out as a unique trigger in their analyses.Based on these considerations, we consider this trigger to be astrophysical and include it in the analysis in the remainder of this paper.

Source Properties
We infer the physical properties of the two GW events using a coherent Bayesian analysis following the methodology described in Appendix B of Abbott et al. (2019c).For GW200105, data from LIGO Livingston and Virgo are analyzed, whereas for GW200115, data from both LIGO detectors and Virgo are used.
Owing to the different signal durations, we analyze 32 s of data for the higher-mass event GW200105 and 64 s of data for GW200115.All likelihood evaluations use a low-frequency cutoff of f low = 20 Hz, except for LIGO Livingston for GW200115, where f low = 25 Hz avoids excess noise localized at low frequencies, as discussed in Section 2. The power spectral density used in the likelihood calculations is the median estimate calculated with BayesLine (Cornish & Littenberg 2015).
The parallel Bilby (PBILBY) inference library, together with the DYNESTY nested sampling software (Ashton et al. 2019;Romero-Shaw et al. 2020b;Smith et al. 2020;Speagle 2020), is the primary tool used to sample the posterior distribution of the sources' parameters and perform hypothesis testing.In addition, we use RIFT (Lange et al. 2018) for the most computationally expensive analyses and LALINFERENCE (Veitch et al. 2015) for verification.
We base our main analyses of GW200105 and GW200115 on BBH waveform models that include the effects of spininduced orbital precession and higher-order multipole GW moments but do not include tidal effects on the secondary.Specifically, we use two signal models: IMRPhenomXPHM (Phenom PHM; Pratten et al. 2020a) from the phenomenological family and SEOBNRv4PHM (EOBNR PHM; Ossokine et al. 2020) from the effective one-body numerical relativity family.The acronym PHM stands for Precessing Higher-order multipole Moments.Henceforth, we will use the shortened names for the waveform models.
The secondary objects are probably NSs based on mass estimates, as discussed in detail in Section 5.As in earlier GW analyses (Abbott et al. 2017a(Abbott et al. , 2020b)), we proceed with two different priors on the secondary's spin magnitude: a low-spin prior, χ 2 0.05, which captures the maximum spin observed in Galactic BNSs that will merge a Hubble time (Burgay et al. 2003), and a high-spin prior, χ 2 0.99, which is agnostic about the nature of the compact object.The two priors allow us to investigate whether the astrophysically relevant subcase of low NS spin leads to differences in the parameter estimation for the binaries.All other priors are set as in previous analyses (e.g., Abbott et al. 2021b).Throughout, we assume a standard flat ΛCDM cosmology with Hubble constant H 0 = 67.9km s −1 Mpc −1 and matter density parameter Ω m = 0.3065 (Ade et al. 2016).
For each spin prior, we run our main analyses with higherorder multipole moments and precession for both waveform families, EOBNR PHM and Phenom PHM.The EOBNR PHM model is used in combination with RIFT and the Phenom PHM model with PBILBY.The parameter estimation results for the individual precessing waveform models yield results in very good agreement; the median values typically differ by 1/10 of the width of the 90% credible interval.Nevertheless, in order to alleviate potential biases due to different samplers or waveform models, we combine an equal number of samples of each into one data set for each spin prior (Abbott et al. 2016c(Abbott et al. , 2020c;;Ashton & Khan 2020) and denote these as Combined PHM.The quoted parameter estimates in the following sections are the Combined PHM high-spin prior analyses.In the figures, we emphasize the high-spin prior results.The values of the most important parameters of the binaries are summarized in Note.We report the median values with 90% credible intervals.Parameter estimates are obtained using the Combined PHM samples.13 The Astrophysical Journal Letters, 915:L5 (24pp), 2021 July 1 Table 2, and we will present the details in the following sections.

Masses
Figure 4 shows the posterior distribution for the component masses of the two binaries.Defining the mass parameters such that the heavier mass is the primary object, i.e., m 1 > m 2 , our analysis shows that GW200105 is a binary with a mass ratio of = = -+ q m m 0.22 The primary components of GW200105 and GW200115 are identified as BHs from their mass measurements.For GW200115, we find that the probability of the primary falling in the lower mass gap (3 M e  m 1  5 M e ; Bailyn et al. 1998;Orosz et al. 1998) is 30% (27%) for a high-spin (low-spin) prior.For context, Figure 4 also includes two potential NSBH candidates discovered previously; GW190814 (Abbott et al. 2020c) is a high-S/N event with well-measured masses that has a significantly more massive primary and a distinctly more massive secondary than either GW200105 or GW200115, and the marginal candidate GW190426_152155 (Abbott et al. 2021b) has (if of astrophysical origin) m 1 -m 2 contours that overlap those of GW200115.The masses of GW190426_152155 are less constrained than those of GW200115 due to its smaller S/N.To highlight how the secondary masses of GW200105 and GW200115 compare to the maximum NS mass, we also show two estimates of the maximum NS mass based on analyses of nonrotating (Landry et al. 2020) and Galactic (Farr & Chatziioannou 2020) NSs.
The secondary masses are consistent with the maximum NS mass, which we quantify in Section 5.2.

Sky Location, Distance, and Inclination
We localize GW200105's source to a sky area of 7200 deg 2 (90% credible region).The large sky area arises due to the absence of data from LIGO Hanford.The luminosity distance of the source is found to be = .For the second event, GW200115, we localize its source to be within 600 deg 2 .It is better localized than GW200105 by an order of magnitude, since GW200115 was observed with three detectors.We find the luminosity distance of the source to be = .The luminosity distance is degenerate with the inclination angle θ JN between the line of sight and the binaries' total angular momentum vector (Cutler & Flanagan 1994;Nissanke et al. 2010).Inclination θ JN = 0 indicates that the angular momentum vector points toward Earth.The posterior distribution of the inclination angle is bimodal and strongly correlated with luminosity distance, as shown in Figure 5.The inclination measurement for GW200105 equally favors orbits that are either oriented toward or away from the line of sight.In contrast, GW200115 shows modest preference for an orientation θ JN π/2.

Spins
The angular momentum vector S i of each compact object is related to its dimensionless spin vector 22 , which is consistent with zero.For GW200115, the spin magnitude is not as tightly constrained, c = -+ 0.33 1 0.29 0.48 , but is also consistent with zero.The spin of the secondary for both events is unconstrained.
One of the best-constrained spin parameters is the effective inspiral spin parameter χ eff (Damour 2001;Racine 2008;Santamaría et al. 2010;Ajith et al. 2014;Vitale et al. 2017b).It encodes information about the binaries' spin components parallel to the orbital angular momentum, , where L is the unit vector along the orbital angular momentum.For GW200105, c = - -+ 0.01 eff 0.15 0.11 , and we find the effective inspiral spin parameter to be strongly peaked about zero, with roughly equal support for being either positive or negative.For GW200115, we find modest support for negative effective inspiral spin: c = - .We find c = - 24 and a probability of 88% that χ 1,z < 0.
The joint posterior probability of the dimensionless spin angular momentum magnitude and tilt angle for both components of both events is shown in Figure 6.The tilt angle with respect to the orbital angular momentum is defined as ( ˆˆ) c L arccos .i .Deviations from uniform shading indicate a spin orientation measurement.The spin orientation of the primary of GW200105 is unconstrained, whereas the orientation of GW200115 shows support for negatively aligned primary spin.
Orbital precession is caused by a spin component in the orbital plane of a binary (Apostolatos et al. 1994), which we parameterize using the effective precession spin parameter 0 χ p 1 (Schmidt et al. 2015).We infer c = -+ 0.09 p 0.07 0.14 for GW200105 and c = -+ 0.21 p 0.17 0.30 for GW200115.To assess the significance of a measurement of precession, we  .57 .For both events and diagnostics, this indicates inconclusive evidence of precession.This result is expected given the S/Ns and inferred inclination angles of the binaries (Vitale et al. 2014;Green et al. 2020;Pratten et al. 2020c).
Low values of the primary mass of GW200115 (m 1  5 M e ) are strongly correlated with negative values of the primary parallel spin component χ 1,z , as shown in Figure 7.The astrophysical implications of the mass and spin correlation are discussed in Section 6. Figure 7 also shows the in-plane spin component χ ⊥ , which is peaked about zero.The lack of conclusive evidence for spin precession in GW200115 is consistent with the measurement of χ ⊥ .Apparent differences between the probability density of the primary spin in Figure 6 and the posteriors of χ 1⊥ -χ 1,z in Figure 7 arise from different choices in visualizing the spin orientation posteriors.

Remnant Properties
Under the hypothesis of NSBH coalescence for the two events, estimates for the final mass and spin of the remnant BH can be made using the models of Zappa et al. (2019) 6. Two-dimensional posterior probability for the spin-tilt angle and spin magnitude for the primary objects (left hemispheres) and secondary objects (right hemispheres) for both events.Spin-tilt angles of 0°(180°) correspond to spins aligned (antialigned) with the orbital angular momentum.The color indicates the posterior probability per pixel of the high-spin prior analysis.For comparison with the low-spin analysis, the solid (dashed) lines indicate the 90% credible regions of the high-spin (low-spin) prior analyses.The tiles are constructed linearly in spin magnitude and cosine of the tilt angles such that each tile contains an identical prior probability.The probabilities are marginalized over the azimuthal angles.
Figure 7. Properties of the primary component of GW200115.The corner plot shows the one-dimensional (diagonal) and two-dimensional (off-diagonal) marginal posterior distributions for the primary's mass and perpendicular and parallel spin components.The shading indicates the posterior probability of the high-spin prior analysis.The solid (dashed) lines indicate the 50% and 90% credible regions of the high-spin (low-spin) prior analyses.The vertical lines indicate the 90% credible intervals for the analyses with high-spin (solid lines) and low-spin (dashed lines) prior.Meidam et al. 2018;Abbott et al. 2019a) show that GW200105 and GW200115 have too low an S/N to allow for tighter constraints than those already presented in Abbott et al. (2021c).Within their measurement uncertainties, our results do not show statistically significant evidence for deviations from the prediction of GR.
To quantify the evidence for higher-order GW multipole moments, we calculate the orthogonal optimal S/N, r lm , for the subdominant multipole moments (Abbott et al. 2020c(Abbott et al. , 2020d;;Mills & Fairhurst 2021).We find (ℓ, |m|) = (3, 3) to be the loudest subdominant multipole moment, as expected for binaries with asymmetric masses.Using the Phenom HM waveform model, we infer r = spin prior.In Gaussian noise, the median of r 33 is approximately χ-distributed with two degrees of freedom, and values greater than 2.1 indicate significant higher-order multipole content.The measured r 33 is therefore consistent with Gaussian noise, as expected for the majority of NSBHs at these S/Ns, except for those viewed close to edge-on.

Waveform Systematics
Our primary results are obtained using precessing BBH models with higher-order multipole moments, Phenom PHM and EOBNR PHM.We now justify this choice by investigating potential systematic uncertainties due to our waveform choice.
First, we investigate the agreement between independent waveform models that incorporate identical physics.Figure 8 shows the two-dimensional m 2 -χ eff posteriors for both events obtained using a variety of NSBH and aligned-spin BBH models.Because some NSBH models only cover χ 1 < 0.5, we restrict the prior range of all models to χ 1 < 0.5 for consistency.
The main panels of Figure 8 are dominated by a correlation of the effective inspiral spin parameter χ eff with the secondary mass m 2 (Cutler & Flanagan 1994;Ng et al. 2018a).Both NSBH models (Phenom NSBH and EOBNR NSBH) give consistent results with each other, as do both BBH models (Phenom and EOBNR), both with and without higher-order multipole moments, with the most notable difference being that EOBNR HM yields tighter posteriors than Phenom HM.This demonstrates that waveform models including the same physics give comparable results, but more studies are warranted to improve the understanding of the BBH waveform models in the NSBH region of parameter space.While not shown in Figure 8, we also find good agreement between the primary precessing BBH waveform models; see Section 4.
Second, comparing the NSBH models with the BBH models without higher-order multipole moments (Phenom and EOBNR), the NSBH models recover similar posterior contours in the m 2 -χ eff plane.This is expected given the asymmetric mass ratio and low S/N of these NSBH observations; see, e.g., Huang et al. (2021) for a demonstration that higher S/Ns would be needed to see notable systematic effects.We observe differences at the extreme ends of the m 2 -χ eff contours (i.e., at the smallest and largest values of m 2 ).The construction of the NSBH waveform models used here did not rely on numerical relativity results at mass ratios q  1/8, nor did they include simulations with χ 1z < 0 or NS masses m 2 > 1.4 M e (Matas et al. 2020;Thompson et al. 2020).Therefore, some differences should be expected, especially for large m 2 in GW200105.Furthermore, for GW200105, the tails of the m 2 -χ eff distribution for Phenom NSBH and EOB NSBH at high m 2 are also impacted by the inability of the data to constrain the tidal deformability.Hence, the posterior samples include combinations of high m 2 with large Λ 2 , despite such combinations being unphysical.This effect is not apparent for GW200115 because of its smaller secondary mass.The isolated islands of probability in the extreme tails of the distributions are due to sampling noise.
Last, when adding the extra physical content of higher-order multipole moments in BBH models (through Phenom HM and EOBNR HM), the extreme ends of the m 2 -χ eff contours are excluded, while the bulk of the distributions are consistent with the posteriors obtained with the NSBH models.In summary, these comparisons indicate that (i) waveform models including the same physics give comparable results; (ii) going from NSBH models to comparable BBH models changes the results only marginally, i.e., any effects of tides are small; and (iii) inclusion of higher-order multipole moments changes the posterior contours more substantially than inclusion of tides.We conclude that the inclusion of precession and higher-order multipole moments afforded by the BBH waveform models is more important than the impact of tides in the NSBH models.

Nature of the Secondary Components
In this section, we describe the investigations to establish the nature of the secondary objects.In Section 5.1, we look for imprints of tidal deformations of the secondaries and conclude that the masses, spins, distances, and S/Ns of the detections make definitive identifications of NSs unlikely in both GW and EM measurements.However, in Section 5.2, we show that the posterior distributions of the secondary masses agree with those of known NSs.

Tidal Deformability and Tidal Disruption
The tidal deformability of NSs is imprinted in the GW signal and investigated using the parameter estimation techniques of Section 4. In contrast, BHs have zero tidal deformability (Binnington & Poisson 2009;Damour & Nagar 2009;Chia 2020;Charalambous et al. 2021;Le Tiec & Casals 2021).We infer the tidal deformability Λ 2 of the NSs in GW200105 and GW200115 using the NSBH waveform models that include tides.We find that the tidal deformabilities are uninformative relative to a uniform prior in Λ 2 ä [0, 5000].This measurement cannot establish the presence of NSs, which is expected given the mass ratios and S/N of the detections (Foucart et al. 2013;Kumar et al. 2017;Fasano et al. 2020;Huang et al. 2020).
Toward the end of the inspiral, the BH may tidally disrupt the NS and form an accretion disk (Pannarale 2013;Foucart et al. 2018).This is hypothesized to drive a relativistic jet (Pannarale et al. 2011;Paschalidis et al. 2015).Given the mass ratios for both events and the aligned spins χ 1,z of their primaries (near zero for GW200105, probably negative for GW200115), we do not expect tidal disruption to occur, which would require more equal masses or more positive χ 1,z (Rantsiou et al. 2008;Shibata & Taniguchi 2008;Etienne et al. 2009;Foucart et al. 2011Foucart et al. , 2018;;Kyutoku et al. 2011;Pannarale 2013).
To quantitatively confirm this expectation, we use the spectral representation of equations of state from Abbott et al. (2018a), which uses an SLY low-density crust model (Douchin & Haensel 2001) and parameterizes the adiabatic index into a polynomial in the logarithm of the pressure (Lindblom 2010;Lindblom & Indik 2012, 2014).Following Stachie et al.
(2021), we marginalize the parameter estimation samples from the NSBH analyses over these equations of state.For a fixed equation of state, we compute the maximum Tolman-Oppenheimer-Volkov (TOV) mass, allowing us to infer the nature (NS or BH) of the lighter binary compact object, as well as its radius R 2 , compactness C 2 = Gm 2 /(R 2 c 2 ), and baryon mass.Based on these, we define a total ejecta mass m ej (Fernández et al. 2019) as the sum of dynamical ejecta (Krüger & Foucart 2020) and 15% of the mass of disk winds (Foucart et al. 2018).For both events, we find that m ej < 10 −6 M e for 99% of the samples.
The absence of ejecta is compatible with the lack of observed EM counterparts.However, given the large distances of the mergers (;300 Mpc) and the large uncertainties of their sky localization, EM emission would have been difficult to detect and associate with these GW events.
Estimating the impact of nonlinear p-g tidal coupling (Abbott et al. 2019b), we find that it would produce a relative frequency-domain phase shift for GW200105 (GW200115) of approximately 134 (38) times smaller than the equivalent phase shift for GW170817.This strongly reduced effect is caused by the larger chirp masses, more asymmetric mass ratios, and the presence of only a single NS.Since p-g effects were not detected for GW170817 (Abbott et al. 2019b), they will be unobservable within the new NSBH systems.

Consistency of Component Masses with the NS Maximum Mass
Even without definite identification of matter signatures in the signals, we can compare the observed m 2 for GW200105 and GW200115 with the maximum NS mass, M max .The existence of massive pulsars (Antoniadis et al. 2013;Cromartie et al. 2019;Fonseca et al. 2021)  (e.g., Shibata et al. 2019;Abbott et al. 2020b;Nathanail et al. 2021).However, rapid rotation could support a larger M max .Given the considerable uncertainty in M max , we examine three different scenarios.Following Essick & Landry (2020), we compute for each scenario the probability ( ) < p m M 2 m a x that the secondary mass is below the maximum NS mass by marginalizing over the uncertainty in M max and the uncertainty of our m 2 measurement.
Supposing that the secondaries are slowly spinning, we consider in the first row of Table 3 an estimate of M max,TOV from a nonparametric astrophysical inference of the equation of state (Landry et al. 2020), which predicts = M max,TOV and is shown in Figure 4. We then relax the lowspin assumption, estimating in the second row the maximum rotationally supported mass ( ) c M max 2 and the breakup spin c max with the universal relations from Breu & Rezzolla (2016).In this scenario, we require and c c  2 m a x for consistency with an NS.Finally, in the third row, we consider a parametric fit to the entire distribution of observed Galactic NSs, including rapidly rotating pulsars (Farr & Chatziioannou 2020), which predicts and is shown in Figure 4.This scenario accounts for the possibility that the maximum mass in the NS population is limited by the astrophysical processes that form compact binaries.Assuming low spin (first row), we find probabilities of 96% and 98% that the secondaries in GW200105 and GW200115, respectively, are consistent with an NS assuming a uniform prior in m 2 .The Note.The values shown use a flat prior in m 2 ; alternative, astrophysically motivated mass priors can cause the estimates to vary by up to 11% across our chosen models.
possibility of large secondary spin reduces these probabilities by up to 3% (second and third rows).So far, this analysis has assumed priors that are uniform in component masses.However, there is considerable uncertainty in the astrophysical mass priors of such systems, and different prior assumptions can affect the component mass posteriors for detections with moderate S/N.To illustrate the impact of population assumptions, we consider three These different mass priors change the numbers in Table 3 by at most 11%, with the smallest values for GW200105 and GW200115 being 89% and 87%, respectively.The decrease is due to the three priors assigning more probability density to equal-mass systems, thus favoring larger m 2 .Thus, the secondaries of both systems are consistent with NSs based on our assumptions about the equation of state and mass distribution.
However, consistency with the maximum NS mass does not exclude the possibility that the secondaries could be BHs or exotic compact objects, if such objects also exist within the NS mass range.For instance, models of primordial BHs predict a peak in the primordial BH mass function at ∼1 M e (Carr et al. 2021).These models also predict that primordial BHs may form coalescing binaries at mass ratios comparable to those reported here.

Astrophysical Implications
The first confident observations of NSBH binaries enable us to study this novel type of astrophysical system in entirely new ways.We pursue three different avenues in this section.First, we infer the merger rate of NSBH binaries in the local universe.
We then place the inferred source properties and merger rate in the context of models of NSBH formation channels and previous EM and GW observations of BHs and NSs.Finally, we investigate to what extent the events reported here can serve to measure the Hubble constant and whether lensing of GWs may have played a role in the observations.

Merger Rate Density
We infer the NSBH merger rate density with our observations using two different approaches.
In the first approach, we consider only GW200105 and GW200115.Following the method of Kim et al. (2003) as previously used in, e.g., Abbott et al. (2016dAbbott et al. ( , 2020c)), we calculate an event-based merger rate assuming one Poisson-distributed count each from GW200105-and GW200115-like populations.We semianalytically calculate the search sensitivities across O1, O2, and the first 9 months of O3 to NSBH populations corresponding to the mass and spin posteriors for the two events.We then calibrate these sensitivities to the results of GSTLAL (Cannon et al. 2012;Privitera et al. 2014;Messick et al. 2017;Sachdev et al. 2019;Hanna et al. 2020) using a broad NSBH-like population and an FAR threshold of 1/(2.8 yr).This FAR threshold is chosen to include only GW200105 and GW200115 while excluding lowsignificance triggers like GW190426_152155, though a more conservative threshold of 1/(100 yr), as used in, e.g., Abbott et al. (2016dAbbott et al. ( , 2020c)), changes the estimated sensitivities by less than 15%.Applying the Poisson Jeffreys prior proportional to - 1 , plotted in green in Figure 9.The second approach to calculating a merger rate takes into account not only GW200105 and GW200115 but also less significant search triggers with masses consistent with the typical range associated with NSBH binaries.Specifically, we consider triggers across O1, O2, and the first 9 months of O3 with associated component masses m 1 ä [2.5,40] M e and m 2 ä [1, 3] M e , i.e., broader ranges than the component mass posteriors for GW200105 and GW200115 obtained in Section 4. The cutoff of m 2 3 M e is chosen as a robust upper limit on the maximum NS mass (Rhoades & Ruffini 1974;Kalogera & Baym 1996).The population is defined to be uniformly distributed in comoving volume, uniform in log component masses in the given ranges, with aligned spins with spin magnitudes distributed uniformly in [0, 0.95].We use as input GSTLAL search triggers above an S/N threshold such that the number of noise triggers exceeds the number of astrophysical signals by a factor of ∼100.Following Farr et al. (2015) and Kapadia et al. (2020), we find the resulting joint likelihood on Poisson parameters for signals of each astrophysical category (Λ BNS , Λ NSBH , Λ BBH ) and for terrestrial noise triggers Λ background .Here the BBH and BNS categories are defined to include triggers where both component masses fall within 5-100 or 1-2.5 M e , respectively.We apply the Jeffreys prior and recover a merger rate density NSBH , where 〈VT〉 NSBH is the population-averaged sensitive time-volume estimated using the GSTLAL pipeline and the NSBH population defined above., the black line in Figure 9.While this rate is higher than our event-based rate, it considers a wider population that includes additional triggers; for example, the component masses of GW190814 (although the nature of its secondary is uncertain) fall within this population.The calculation is also based on the component mass values of GSTLAL search triggers, adjusted with an analytical model (Fong 2018) of the response of the template bank to signals with a given distribution of true masses.This procedure is expected to be less precise than Bayesian parameter estimation, which is impractical for the large numbers of triggers involved.
The merger rate density measured here is consistent with the upper bound of 610 Gpc −3 yr −1 derived from the absence of NSBH detections in O1 and O2 (Abbott et al. 2019c).Revised merger rate estimates for all CBC sources will be provided in a future analysis of the full O3 data set.

System Origins
To understand the origin of GW200105 and GW200115, we compare their observed masses and spins with theoretical predictions.Population synthesis studies modeling the various formation channels of merging compact object binaries distinguish between NSs and BHs with a simple mass cut, typically between 2 and 3 M e .This is consistent with the secondary masses of both events being classified as NSs, so for the purposes of discussing formation channels for these events and rates, we will discuss GW200105 and GW200115 in the context of them being NSBHs.

Formation Channels
Formation channels for NSBHs can be broadly categorized as either isolated binary evolution or one of several dynamical formation channels (e.g., globular clusters or nuclear star clusters).Since isolated binaries form in young star clusters and can be influenced by dynamical interactions before the cluster dissolves and the binary effectively becomes isolated, rates from young star clusters naturally encompass rates from isolated binaries.Predicted rates of NSBH mergers in the local universe vary by orders of magnitude across the various formation channels.
Models of the canonical isolated binary evolution channelin which stellar progenitors evolve together, shedding orbital angular momentum through phases of stable and/or unstable mass transfer prior to compact object formation-predict NSBH merger rate densities around 0.1-800 Gpc −3 yr −1 (Belczynski et al. 2002;Sipior & Sigurdsson 2002;Belczynski et al. 2006Belczynski et al. , 2016;;Dominik et al. 2015;Eldridge et al. 2017;Kruckow et al. 2018;Mapelli & Giacobbo 2018;Neijssel et al. 2019;Drozda et al. 2020;Zevin et al. 2020;Broekgaarden et al. 2021).The high uncertainty is driven by the lack of observed NSBHs and the wide range of model assumptions.Merger rates are sensitive to the treatment of common envelopes, which may be a necessary evolutionary phase for producing compact binaries in tight orbits capable of merging in a Hubble time (Ivanova et al. 2013).They are also sensitive to prescriptions for supernova kick magnitudes.While moderate kicks can produce eccentric orbits that merge on short timescales, high supernova kicks may disrupt the progenitor binaries and suppress the merger rate (Belczynski et al. 2002;Giacobbo & Mapelli 2020;Tang et al. 2020).
Models of dynamical formation channels in denser environments typically predict much lower merger rates.For instance, in globular clusters and nuclear star clusters, BHs segregate and dominate the core, where the bulk of dynamical interactions occur (Portegies Zwart & McMillan 2000;Morscher et al. 2015), so that encounters between NSs and BHs are relatively rare, with NSBH merger rate densities on the order of 10 −2 Gpc −3 yr −1 (Clausen et al. 2013;Arca Sedda 2020;Ye et al. 2020).In disks of active galactic nuclei, the presence of gas could possibly increase the NSBH merger rate density up to 300 Gpc −3 yr −1 (McKernan et al. 2020).
NSBHs may also merge via hierarchical triple interactions, where inner NSBH binaries are driven to high eccentricity by massive tertiary companions and merge on rapid timescales (Antonini et al. 2017;Silsbee & Tremaine 2017).However, the predicted merger rates are negligible unless supernova kicks are assumed to be zero (Fragione & Loeb 2019).
It is likely that a combination of the above channels contribute to the astrophysical NSBH merger rate.However, the isolated binary evolution, young star cluster, and active galactic nuclei channels are capable of individually accounting for the NSBH merger rate estimated here.

Masses
While there are no observed NSBHs in the Milky Way, we can place the component masses of GW200105 and GW200115 in the context of the observed population of BH and NS masses, as well as the predicted populations of NSBHs.Observations suggest that the mass distribution of the Galactic population of NSs peaks around 1.33 M e , with a secondary peak around 1.9 M e (Antoniadis et al. 2016;Alsing et al. 2018).The secondary mass observed in GW200115 and marginal event GW190426_152155 are consistent with the population peaking at 1.33 M e , while the secondary observed in GW200105 (;1.9 M e ) and the primary component from BNS merger GW190425 (m 1 = 1.60-1.87M e ; Abbott et al. 2020b) are consistent with the high-mass population.However, a rigorous association of the events with different components of the NS population would require a thorough population analysis.Radio observations of BNS systems do not find such massive NSs, leading to speculation as to the origin of GW190425 (Romero-Shaw et al. 2020a;Safarzadeh et al. 2020;Galaudage et al. 2021;Mandel et al. 2021).Stellar metallicities in the Milky Way are not representative of all populations of GW sources (O'Shaughnessy et al. 2008(O'Shaughnessy et al. , 2010;;Belczynski et al. 2010;Eldridge et al. 2017;Neijssel et al. 2019).
The BH masses observed in GW200105 and GW200115 ( - + 8.9 1.5 1.2 and , respectively) are in line with predictions from population synthesis models for NSBH mergers from isolated binary evolution and young star clusters.In NSBHs, current binary evolution models do not predict BH masses above ;10 M e (Eldridge et al. 2017;Kruckow et al. 2018;Mapelli & Giacobbo 2018;Broekgaarden et al. 2021), while Population III NSBHs (Kinugawa et al. 2017) and dynamical interactions in low-metallicity young star clusters allow for higher BH masses (Ziosi et al. 2014;Rastello et al. 2020).
Electromagnetic observations of X-ray binaries have not uncovered BHs between 3 and 5 M e , leading to speculation about a mass gap (Özel et al. 2010;Farr et al. 2011;Kreidberg et al. 2012;Miller & Miller 2014).Analysis of GWTC-2 has also found evidence for a gap or dip in the BH mass spectrum between ∼2.6 and 4 M e (Fishbach et al. 2020).For GW200115, we find nonnegligible support for the primary lying in this mass gap, with p(3 M e < m 1 < 5 M e ) = 30% (27%) under the high-spin (low-spin) priors.This low-mass region is correlated with negative values of the parallel component of the primary spin; see Section 6.2.3 below.
In summary, the masses inferred for GW200105 and GW200115 are consistent with expectations for NSBHs; their primary masses are in agreement with predictions for BH masses in population synthesis models of the dominant formation scenarios.Meanwhile, their secondary masses are compatible with the observed population of Galactic NSs, as well as the masses inferred from GW observations of BNS mergers.
While the secondary spins of both events reported here are poorly constrained due to the unequal masses, the primary spins of GW200105 and GW200115 can be placed in the context of predictions for BH spins from models of stellar and dynamical evolution and EM observations of NSBH progenitors.As can be seen in Figures 6 and 7, the primary spins of GW200105 and GW200115 are consistent with zero ( -+ 0.08 0.08 0.22 and -+ 0.33 0.29 0.48 , respectively), but moderate values of spin are not ruled out.The primary in GW200115 may even have relatively high spin, with a 90% upper limit of 0.72.Several studies of the observed population of high-mass X-ray binaries (Liu et al. 2008;Gou et al. 2009Gou et al. , 2014;;Zhao et al. 2021) find that the BHs have large spins (Valsecchi et al. 2010;Wong et al. 2012;Qin et al. 2019;Zhao et al. 2021).Given the short lifetime of the secondary, mass transfer is argued to be insufficient to generate BHs with such high spins, implying that the BHs were born with high spins.Belczynski et al. (2011) found that one such high-mass X-ray binary, Cygnus X-1, is expected to form an NSBH with a BH that carries near-maximal spin, although it would not merge within a Hubble time.However, following revised estimates of the component masses of Cygnus X-1 (Miller-Jones et al. 2021), Neijssel et al. (2021) found that it will most likely form a BBH.Meanwhile, analyses of GWTC-1 and GWTC-2 have found evidence for BH spin (Abbott et al. 2016f, 2021dVitale et al. 2017a;Chatziioannou et al. 2019;Kimball et al. 2020), though they do not determine whether those BHs may have been formed with that spin.Altogether, these EM and GW observations of compact binaries and their progenitors suggest a range of BH natal spins in NSBH binaries.
Along with their magnitudes, the alignments of component spins with the overall binary orbital angular momentum are of astrophysical interest.In particular, we find evidence that the primary BH spin in GW200115 is negatively aligned with respect to the orbital angular momentum axis, with p(χ 1,z < 0) = 88% (87%) under the high-spin (low-spin) prior and the more negative values of χ 1,z correlated with smaller m 1 .This negative alignment is consistent with dynamical formation channels, which typically form binaries with random spin orientations (Rodriguez et al. 2016), but the predicted rates from these channels, discussed in Section 6.2, are small.Binaries born in isolation are expected to form with only small misalignments (30°; Kalogera 2000), though they may become misaligned by supernova kicks at compact object formation (Rodriguez et al. 2016;Gerosa et al. 2018;Wysocki et al. 2018) and possibly during subsequent evolution via mass transfer (Stegmann & Antonini 2021).Meanwhile, NSBH progenitor binaries originating in young star clusters can be perturbed via close dynamical encounters before being ejected into the field.Therefore, a misaligned spin in the primary of GW200115 would not necessarily rule out any of the plausible NSBH formation channels.

Cosmology and Lensing
Gravitational wave sources are standard sirens, providing a direct measurement of their luminosity distance (Schutz 1986;Holz & Hughes 2005), and they can be used to measure the Hubble constant (Abbott et al. 2017c(Abbott et al. , 2021a;;Fishbach et al. 2019).Due to the lack of a confirmed EM counterpart and large numbers of galaxies inside the localization volumes of each of the two events, we do not obtain any informative bounds on H 0 from these observations.The detections of GW200105 and GW200115 are separated by only ∼10 days.One explanation for the small time delay could be that the two events are created by gravitational lensing by a galaxy (Haris et al. 2018).Gravitational lensing is unlikely even a priori (Smith et al. 2017;Ng et al. 2018b), and the nonoverlapping mass posteriors (Figure 4) further exclude it as a possible explanation (Haris et al. 2018).While GW200115 and GW190426_152155 exhibit agreement in their source mass posteriors, their sky localization areas do not overlap, and their detector-frame (redshifted) chirp masses show only marginal overlap (Abbott et al. 2021b), ruling out lensing as a possible explanation.

Conclusions
During its third observing run, the LIGO-Virgo GW detector network observed GW200105 and GW200115, two GW events consistent with NSBH coalescences.Event GW200105 is effectively a single-detector event observed in LIGO Livingston with an S/N of 13.9.It clearly stands apart from all recorded noise transients, but its statistical confidence is difficult to establish.Event GW200115 was observed in coincidence by the network with an S/N of 11.6 and FAR of <1/(1 × 10 5 yr).
The source component masses of GW200105 and GW200115 make it likely that these events originated from NSBH coalescences.Their primary masses are found to be = -+ m 8.9 , respectively, are consistent with the observed NS mass distribution in the Milky Way, as well as population synthesis predictions for secondary masses in merging NSBHs.
We find no evidence of measurable tides or tidal disruption for either of the two signals, and no EM counterparts to either detection have been identified.As such, there is no direct evidence that the secondaries are NSs, and our observations are consistent with either event being a BBH merger.However, the absence of tidal measurements and EM counterparts is to be expected given the properties and distances of the two events.Moreover, the comparisons of the secondary masses to the maximum allowed NS mass yield a probability ( )  p m M 2 m a x of 89%-96% and 87%-98% for the secondaries in GW200105 and GW200115, respectively, being compatible with NSs (see Section 5.2).The effective inspiral spin parameter of GW200105 is strongly peaked around zero: c = - 24 , and we find support for negatively aligned primary spin (χ 1,z < 0) at 88% probability.More negative values of χ 1,z in GW200115 are correlated with particularly small primary masses reaching into the lower mass gap.We find p(3 M e < m 1 < 5 M e ) = 30% (27%) under the high-spin (low-spin) parameter estimation priors.We find no conclusive evidence for spin-induced orbital precession in either system.
We estimate the merger rate density of NSBH binaries with two approaches.Assuming that GW200105 and GW200115 are representative of the entire NSBH population, we find 1 .These are the first direct measurements of the NSBH merger rate.Both estimates are broadly consistent with the rate predicted from NSBH formation in isolated binaries or via young star clusters.Formation channels in dense star clusters (globular or nuclear) and triples predict lower rates than those inferred from the two events and are unlikely to be the dominant NSBH formation channels.
The observations of GW200105 and GW200115 are consistent with predictions for merging NSBHs and observations of BHs and NSs in the Milky Way.Given their significantly unequal component masses, future observations of NSBH systems will provide new opportunities to study matter under extreme conditions, including tidal disruption, and search for potential deviations from GR.
Gravitational Physics (Hypatia), and the Barcelona Supercomputing Center-Centro Nacional de Supercomputación.
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.
may have been an NSBH merger.Although GW190814ʼs secondary mass of mass supported by slowly spinning NSs

Figure 1 .
Figure 1.Time-frequency representations of the data containing GW200105 (left column) and GW200115 (right column).Times are shown relative to the signals' merger times, 2020 January 5 at 16:24:26 UTC (left) and 2020 January 15 at 04:23:10 UTC (right).The amplitude scale of each time-frequency tile is normalized by the respective detector's noise amplitude spectral density.The LIGO Livingston data of GW200105 show a track of excess power with increasing frequency.In the other panels, no similar tracks are visible, as the S/N in each of the detectors is lower and (for GW200115) the signal is longer.For GW200105, the LIGO Livingston data are shown after glitch subtraction.For GW200115, light-scattering noise is visible in LIGO Livingston below 25 Hz.

Figure 2 .
Figure 2. Sky localizations for GW200105 (top) and GW200115 (bottom) in terms of R.A. and decl.The thick solid contours show the 90% credible regions from the low-latency sky localization algorithm BAYESTAR (Singer & Price 2016).The shaded patch is the sky map obtained from the preferred high-spin analysis of Section 4, with the 90% credible regions bounded by the thin dotted contours.

Figure 3 .
Figure 3. Colored shading shows the joint S/N-ξ 2 noise probability density function for LIGO Hanford (LHO), LIGO Livingston (LLO), and Virgo, computed from background triggers found by GSTLAL in the region of lighter binary systems during the entire O3.The red star indicates GW200105.At its position, there is no background present, indicating that it stands above all of the recorded background triggers.For comparison, the triggers of GW200115 and the marginal GW190426_152155 are also shown.

Figure 4 .
Figure 4. Component masses of GW200105 (red) and GW200115 (blue), represented by their two-and one-dimensional posterior distributions.Colored shading and solid curves indicate the high-spin prior, whereas dashed curves represent the low-spin prior.The contours in the main panel, as well as the vertical and horizontal lines in the top and right panels, respectively, indicate the 90% credible intervals.Also shown in gray are two possible NSBH events, GW190814 and the marginal candidate GW190426_152155, the latter overlapping GW200115.Lines of constant mass ratio are indicated in dashed gray.The green shaded curves in the right panel represent the one-dimensional probability densities for two estimates of the maximum NS mass, based on analyses of nonrotating NSs (M ; max,TOV Landry et al. 2020, 2021) and Galactic NSs (M ; max,GNS Farr & Chatziioannou 2020).

Figure 5 .
Figure 5. Two-and one-dimensional posterior distributions for distance D L and inclination θ JN .The solid (dashed) lines indicate the high-spin (low-spin) prior analysis, and the shading indicates the posterior probability of the highspin prior analysis.The contours in the main panel and the horizontal lines in the right panel indicate 90% credible intervals.
We use samples obtained by combining those from Phenom NSBH and EOB NSBH.For GW200105, the remnant mass and spin are .We do not investigate any postmerger GW signals.The S/Ns of GW200105 and GW200115 are around a factor of 3 less than that of GW170817, for which there was no evidence of GWs after the merger(Abbott et al. 2017b).In the absence of tidal disruption, the postmerger signals of GW200105 and GW200115 would likely resemble a BH ringdown(Foucart et al. 2013).The GW signal associated with such ringdowns would appear well outside of LIGO's and Virgoʼs sensitive bandwidth given the remnant masses and spins of the systems (Sarin & Lasky 2021).4.5.Tests of General Relativity and Higher-order GW Multipole MomentsResults from parameterized tests of general relativity (GR;Blanchet & Sathyaprakash 1995;Yunes & Pretorius 2009;Mishra et al. 2010;Li et al. 2012aLi et al. , 2012b;;Agathos et al. 2014; ) for GW200115 with the low (high)

Figure 8 .
Figure8.Comparison of two-dimensional m 2 -χ eff posteriors for the two events reported here, using various NSBH and BBH signal models.The vertical dashed lines indicate several mass-ratio references mapped to m 2 for the median estimate of the chirp masses of GW200105 and GW200115.
mass.Studies of GW170817ʼs remnant typically suggest that the maximum mass of a nonrotating NSthe TOV mass-is alternative priors: one based on Salpeter mass distributions, p(m) ∼ m −2.3 (Salpeter 1955), independently for each component; one based on an extrapolation of the BBH mass model BROKEN POWER LAW from Abbott et al. (2021d) down to 0.5 M e for both components; and another based on a similar extrapolation of the POWER LAW + PEAK BBH mass model from the same study.We marginalize over the uncertainties in the latter two models, which are fit to the BBH population from Abbott et al. (2021b), including the outlier event GW190814 with a secondary component mass below 3 M e (Abbott et al. 2020c).
Gpc −3 yr −1 .The PYCBC detection of GW200115, using the same method but with an independent set of injections for calibration(Tiwari 2018), yields a consistent = yr −1 .Combining the likelihoods over  200105 and  200115 according toKim et al. (2003) and applying the Jeffreys prior to the total rate, we find the total event-based NSBH merger rate density =

Figure 9 .
Figure 9. Inferred probability densities for the NSBH merger rate.Green line: rate assuming one count each from a GW200105-and GW200115-like NSBH population.Black line: rate for a broad NSBH population with a low threshold that accounts for marginal triggers.The short vertical lines indicate the 90% credible intervals.

.
For the second event, GW200115, the effective inspiral spin parameter is inferred to be c = - .For GW200115, the spin component parallel to the orbital angular momentum of the primary is c = -

75
Gpc −3 yr −1 .Conversely, assuming a broader range of allowed primary and secondary masses, and considering all triggers in O3, we find =

Table 1
Network S/N Recovered by the Search Pipelines in Low Latency and Later Offline Analysis with Improved Calibration and Refined Data Quality Information

Table 3
Probability that the Secondary Mass Is below the Maximum NS Mass M max for Each Event, Given Different Spin Assumptions and Different Choices for the Maximum NS Mass