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

Event GW200115 was initially identified by all four low-latency matched-filtering pipelines as a possible compact binary coalescence (CBC) candidate in LIGO Hanford and LIGO Livingston: GstLAL (Cannon et al. 2012; Privitera et al. 2014; Messick et al. 2017; Sachdev et al. 2019; Hanna et al. 2020), MBTA Online (Adams et al. 2016; Aubin et al. 2021), PyCBC Live (Usman et al. 2016; Nitz et al. 2017, 2018, 2019a; Dal Canton et al. 2020), and SPIIR (Hooper et al. 2012; Luan et al. 2012; Guo et al. 2018; Chu et al. 2020). As detailed in Abbott et al. (2021b), matched-filtering searches use banks (Owen & Sathyaprakash 1999; Harry et al. 2009; Privitera et al. 2014; Roy et al. 2017, 2019) of modeled gravitational waveforms with the mass range relevant to GW200105 and GW200115 covered by the SEOBNRv4_ROM waveform model (Pürrer 2016; Bohé et al. 2017). The signal in Virgo was not loud enough to further improve the significance of the coincident candidate observed by the LIGO detectors, so GW200115 was reported as a two-detector event in the GraceDB event database named S200115j with a merger time of 2020 January 15 at 04:23:10 UTC.

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 low-latency 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². The GCN Circular reported a probability >99% for the lighter compact object to have a mass below 3 M_⊙, 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, 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, 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⁵ 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⁴ yr), respectively. Based on two detectors, these FARs robustly indicate the confidence of the detection.

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² 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_⊙. The colored region indicates the density of background noise triggers quantified through the S/N–ξ² noise probability density function for the three detectors, where the autocorrelation ξ² provides a consistency test similar to χ² (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.

Zoom In Zoom Out Reset image size

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 single-detector 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.

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₁ > m₂, our analysis shows that GW200105 is a binary with a mass ratio of $q={m}_{2}/{m}_{1}={0.22}_{-0.04}^{+0.08}$, with source component masses ${m}_{1}={8.9}_{-1.5}^{+1.2}$ and ${m}_{2}={1.9}_{-0.2}^{+0.3}\,{M}_{\odot }$. Similarly, GW200115 is a binary with a mass ratio of $q={0.26}_{-0.10}^{+0.35}$, with source component masses ${m}_{1}={5.7}_{-2.1}^{+1.8}$ and ${m}_{2}={1.5}_{-0.3}^{+0.7}\,{M}_{\odot }.$

Zoom In Zoom Out Reset image size

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_⊙ ≲ m₁ ≲ 5 M_⊙; 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₁–m₂ 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.

We localize GW200105's source to a sky area of 7200 deg² (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 ${D}_{{\rm{L}}}={280}_{-110}^{+110}\,\mathrm{Mpc}$. For the second event, GW200115, we localize its source to be within 600 deg². 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 ${D}_{{\rm{L}}}={300}_{-100}^{+150}\,\mathrm{Mpc}$.

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.

Zoom In Zoom Out Reset image size

The angular momentum vector S _i of each compact object is related to its dimensionless spin vector ${{\boldsymbol{\chi }}}_{i}\equiv c{{\boldsymbol{S}}}_{i}/({{Gm}}_{i}^{2})$. Its magnitude χ_i ≡ ∣ χ _i ∣ is bounded by 1. For GW200105, we infer ${\chi }_{1}={0.08}_{-0.08}^{+0.22}$, which is consistent with zero. For GW200115, the spin magnitude is not as tightly constrained, ${\chi }_{1}={0.33}_{-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, ${\chi }_{\mathrm{eff}}=\left(\tfrac{{m}_{1}}{M}{{\boldsymbol{\chi }}}_{1}+\tfrac{{m}_{2}}{M}{{\boldsymbol{\chi }}}_{2}\right)\cdot \hat{L}$, where $\hat{L}$ is the unit vector along the orbital angular momentum.

For GW200105, ${\chi }_{\mathrm{eff}}=-{0.01}_{-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: ${\chi }_{\mathrm{eff}}=-{0.19}_{-0.35}^{+0.23}$. Negative values of χ_eff indicate binaries with at least one spin component negatively aligned with respect to the orbital angular momentum, i.e., ${\chi }_{i,z}\equiv {{\boldsymbol{\chi }}}_{i}\cdot \hat{L}\lt 0$. We find ${\chi }_{1,z}=-{0.19}_{-0.50}^{+0.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 $\arccos \left(\hat{L}.{\hat{\chi }}_{i}\right)$. 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.

Zoom In Zoom Out Reset image size

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 ${\chi }_{{\rm{p}}}={0.09}_{-0.07}^{+0.14}$ for GW200105 and ${\chi }_{{\rm{p}}}={0.21}_{-0.17}^{+0.30}$ for GW200115. To assess the significance of a measurement of precession, we compute a Bayes factor between a precessing and nonprecessing signal model and the precession S/N ρ_p (Fairhurst et al. 2020a, 2020b). For GW200105, we find a log Bayes factor in favor of spin precession of ${\mathrm{log}}_{10}{ \mathcal B }=-0.24$ and precession S/N ${\rho }_{{\rm{p}}}={0.74}_{-0.61}^{+1.35}$. For GW200115, ${\mathrm{log}}_{10}{ \mathcal B }=-0.12$ and ${\rho }_{{\rm{p}}}={0.97}_{-0.79}^{+1.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₁ ≲ 5 M_⊙) 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.

Zoom In Zoom Out Reset image size

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). We use samples obtained by combining those from Phenom NSBH and EOB NSBH. For GW200105, the remnant mass and spin are ${M}_{{\rm{f}}}={10.4}_{-2.0}^{+2.7}$ and ${\chi }_{{\rm{f}}}={0.43}_{-0.03}^{+0.04}$, while for GW200115, ${M}_{{\rm{f}}}={7.8}_{-1.6}^{+1.4}$ and ${\chi }_{{\rm{f}}}={0.38}_{-0.02}^{+0.04}$. 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).

Results from parameterized tests of general relativity (GR; Blanchet & Sathyaprakash 1995; Yunes & Pretorius 2009; Mishra et al. 2010; Li et al. 2012a, 2012b; Agathos et al. 2014; 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, ${\rho }_{{lm}}^{\perp }$, for the subdominant multipole moments (Abbott et al. 2020c, 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 ${\rho }_{33}^{\perp }={1.70}_{-1.09}^{+0.92}$ (${1.70}_{-1.11}^{+0.94}$) for GW200105 and ${\rho }_{33}^{\perp }={0.91}_{-0.66}^{+0.93}$ (${0.86}_{-0.65}^{+0.90}$) for GW200115 with the low (high) spin prior. In Gaussian noise, the median of ${\rho }_{33}^{\perp }$ is approximately χ-distributed with two degrees of freedom, and values greater than 2.1 indicate significant higher-order multipole content. The measured ${\rho }_{33}^{\perp }$ 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.

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₂–χ_eff posteriors for both events obtained using a variety of NSBH and aligned-spin BBH models. Because some NSBH models only cover χ₁ < 0.5, we restrict the prior range of all models to χ₁ < 0.5 for consistency.

Zoom In Zoom Out Reset image size

The main panels of Figure 8 are dominated by a correlation of the effective inspiral spin parameter χ_eff with the secondary mass m₂ (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₂–χ_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₂–χ_eff contours (i.e., at the smallest and largest values of m₂). 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₂ > 1.4 M_⊙ (Matas et al. 2020; Thompson et al. 2020). Therefore, some differences should be expected, especially for large m₂ in GW200105. Furthermore, for GW200105, the tails of the m₂–χ_eff distribution for Phenom NSBH and EOB NSBH at high m₂ are also impacted by the inability of the data to constrain the tidal deformability. Hence, the posterior samples include combinations of high m₂ with large Λ₂, 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₂–χ_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.

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 Λ₂ 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 Λ₂ ∈ [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. 2011, 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₂, compactness C₂ = Gm₂/(R₂ c²), 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⁻⁶ M_⊙ 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.

Even without definite identification of matter signatures in the signals, we can compare the observed m₂ 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) places a lower bound of ${M}_{\max }\gtrsim 2\,{M}_{\odot }$ on the maximum NS mass. Studies of GW170817's remnant typically suggest that the maximum mass of a nonrotating NS—the TOV mass—is ${M}_{\max ,\mathrm{TOV}}\lesssim 2.3\,{M}_{\odot }$ (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}_{2}\lt {M}_{\max })$ that the secondary mass is below the maximum NS mass by marginalizing over the uncertainty in ${M}_{\max }$ and the uncertainty of our m₂ measurement.

Supposing that the secondaries are slowly spinning, we consider in the first row of Table 3 an estimate of ${M}_{\max ,\mathrm{TOV}}$ from a nonparametric astrophysical inference of the equation of state (Landry et al. 2020), which predicts ${M}_{\max ,\mathrm{TOV}}\,={2.22}_{-0.20}^{+0.30}\,{M}_{\odot }$ and is shown in Figure 4. We then relax the low-spin assumption, estimating in the second row the maximum rotationally supported mass ${M}_{\max }({\chi }_{2})$ and the breakup spin ${\chi }_{\max }$ with the universal relations from Breu & Rezzolla (2016). In this scenario, we require ${m}_{2}\leqslant {M}_{\max }({\chi }_{2})$ and ${\chi }_{2}\leqslant {\chi }_{\max }$ 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 ${M}_{\max ,\mathrm{GNS}}\,={2.25}_{-0.31}^{+0.71}\,{M}_{\odot }$ 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₂. The possibility of large secondary spin reduces these probabilities by up to 3% (second and third rows).

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

		$p({m}_{2}\lt {M}_{\max })$
Spin Prior	Choice of ${M}_{\max }$	GW200105	GW200115
∣χ2∣ < 0.05	${M}_{\max ,\mathrm{TOV}}$	96%	98%
∣χ2∣ < 0.99	${M}_{\max }({\chi }_{2})$	94%	95%
∣χ2∣ < 0.99	${M}_{\max ,\mathrm{GNS}}$	93%	96%

Note. The values shown use a flat prior in m₂; alternative, astrophysically motivated mass priors can cause the estimates to vary by up to 11% across our chosen models.

Download table as:  ASCIITypeset image

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 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_⊙ 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_⊙ (Abbott et al. 2020c).

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₂. 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_⊙ (Carr et al. 2021). These models also predict that primordial BHs may form coalescing binaries at mass ratios comparable to those reported here.

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. (2016d, 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 low-significance triggers like GW190426_152155, though a more conservative threshold of 1/(100 yr), as used in, e.g., Abbott et al. (2016d, 2020c), changes the estimated sensitivities by less than 15%. Applying the Poisson Jeffreys prior proportional to ${{ \mathcal R }}_{i}^{-1/2}$, we find ${{ \mathcal R }}_{200105}={16}_{-14}^{+38}$ Gpc⁻³ yr⁻¹ and ${{ \mathcal R }}_{200115}\,={36}_{-30}^{+82}$Gpc⁻³ yr⁻¹. The PyCBC detection of GW200115, using the same method but with an independent set of injections for calibration (Tiwari 2018), yields a consistent ${{ \mathcal R }}_{200115}={40}_{-34}^{+92}$ Gpc⁻³ yr⁻¹. Combining the likelihoods over ${ \mathcal R }$ ₂₀₀₁₀₅ and ${ \mathcal R }$ ₂₀₀₁₁₅ according to Kim et al. (2003) and applying the Jeffreys prior to the total rate, we find the total event-based NSBH merger rate density ${{ \mathcal R }}_{\mathrm{NSBH}}={45}_{-33}^{+75}\,{\mathrm{Gpc}}^{-3}\,{\mathrm{yr}}^{-1}$, plotted in green in Figure 9.

Zoom In Zoom Out Reset image size

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₁ ∈ [2.5, 40] M_⊙ and m₂ ∈ [1, 3] M_⊙, i.e., broader ranges than the component mass posteriors for GW200105 and GW200115 obtained in Section 4. The cutoff of m₂ ≤ 3 M_⊙ 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_⊙, respectively. We apply the Jeffreys prior and recover a merger rate density ${{ \mathcal R }}_{\mathrm{NSBH}}={{\rm{\Lambda }}}_{\mathrm{NSBH}}/\langle {VT}{\rangle }_{\mathrm{NSBH}}$, where 〈VT〉_NSBH is the population-averaged sensitive time–volume estimated using the GstLAL pipeline and the NSBH population defined above. This yields an estimated NSBH merger rate density of ${{ \mathcal R }}_{\mathrm{NSBH}}={130}_{-69}^{+112}\,{\mathrm{Gpc}}^{-3}\,{\mathrm{yr}}^{-1}$, 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⁻³ yr⁻¹ 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.

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_⊙. 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 predicted rates, we will discuss GW200105 and GW200115 in the context of them being NSBHs.

6.2.1. Formation Channels

6.2.2. Masses

6.2.3. Spins

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, 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₀ 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.