Distinguishing binary neutron star from neutron star-black hole mergers with gravitational waves
Hsin-Yu Chen Katerina Chatziioannou












The gravitational-wave signal from the merger of two neutron stars cannot be easily differentiated from 
the signal produced by a comparable-mass mixed binary of a neutron star and a black hole. 
Indeed, both binary types can account for the gravitational-wave signal GW170817 
even if its electromagnetic counterpart emission is taken into account.
We propose a method that requires neither information from the post-inspiral phase of the binary nor an electromagnetic counterpart to identify mixed binaries of neutron stars merging with low-mass black holes 
using gravitational-waves alone. 
This method is based on the fact that certain neutron star properties that can be measured with gravitational-waves
are common or similar for all neutron stars.
For example all neutron stars share the same equation of state and if the latter is hadronic, neutron stars have similar radii.
If a mixed binary is misidentified as a neutron star binary, the inferred neutron star 
properties will be misestimated and appear as outliers in a population of low-mass binaries.
We show that as few as ∼ 5 low-mass events will allow for the identification of the type of one 
event at the 80% confidence level. 
We model the population of low-mass binaries with a hierarchical mixture model
and show that we can constrain the existence of mixed binaries or 
measure their abundance relative to neutron star binaries to 
∼ 0.1 at the 68% credible level with 100 events. 

 



Distinguishing binary neutron star from neutron star-black hole mergers with gravitational waves
    Katerina Chatziioannou
    December 30, 2023
================================================================================================




§ INTRODUCTION



The gravitational-wave (GW) event GW170817 detected by Advanced LIGO <cit.> and Virgo <cit.>
 is consistent with the merger of two neutron stars (BNS) <cit.>. 
 Although the GW data place a lower limit on the compactness 
of the two coalescing bodies, objects more compact than neutron stars (NSs) are not ruled out <cit.>. 
Arriving after the GW signal, the electromagnetic (EM) counterparts GRB 170817A <cit.> 
and kilonova AT 2017gfo, e.g. <cit.>,
imply the presence of at least one NS in the binary. 
However, we still can not exclude the possibility of GW170817 being a merger of a NS and a black hole (NSBH) <cit.>.













X-ray binaries suggest a lack of BHs with mass below 5M_⊙ <cit.>, 
but the origin of this mass gap between BHs and NSs is not fully understood <cit.>. 
Recently, the discovery of a 3.3^+2.8_-0.7M_⊙ unseen companion of the giant star 2MASS J05215658+4359220 
further challenged the existence of the mass gap <cit.>.
Scenarios for the production of low-mass BHs 
include primordial density fluctuations <cit.>, slow supernova explosions <cit.>, 
 mergers of NSs <cit.>, and interactions of dark matter and NSs <cit.>.
 Low-mass binary mergers can potentially help 
study the black hole mass distribution <cit.>, but 
probing the existence of objects in the mass gap is challenging <cit.>.

As already noted in <cit.>, though, constraining the component 
masses in ∼ (0.5-2) M_⊙ does not 
definitively prove the type of the binary. 

For that we also need to detect (or rule out) tidal interactions in the binary with GWs,
 quantified through the NS tidal deformabilities <cit.>.
For binary mergers, the individual tidal parameter of each star is difficult to measure;
instead constraints are placed on a combination of masses and tidal deformabilities, 
Λ̃ <cit.>. For an NSBH that is not particularly loud (signal-to-noise 
ratio less than ≃ 30), the tidal deformability is generally difficult to 
measure <cit.>. 

GW170817 data place a lower limit on Λ̃ subject to the assumption of 
small spins <cit.>; the data are nonetheless consistent with 
a highly spinning BH binary. At the same time, a nonzero Λ̃ only suggests the presence of
one NS, still allowing for the NSBH scenario. 
Further analysis of the EM counterpart remains inconclusive and cannot rule out the NSBH scenario <cit.>.
Similar analyses for near-future detections are subject to the availability and interpretation of an EM counterpart,
while post-merger information <cit.> or evidence for disruption <cit.>
will likely be buried in detector noise.

The misidentification of a low-mass NSBH for a BNS
can have dire consequences for our ability to accurately measure the radius of NSs with GWs. Indeed
a GW analysis of a NSBH assuming it is a BNS underestimates the true radius <cit.>. 
The amount of bias depends on the mass of the BH as the tidal deformability
is a steeply decreasing function of the mass. Misidentifying a ∼ 2 M_⊙ BH
for a NS induces a negligible error, while misidentifying a ∼ 1 M_⊙ BH
can lead to a radius error of multiple km.

We present a method to distinguish between BNSs and low-mass NSBHs using their
GW signals alone. We take advantage of the inferred radius bias that is incurred for NSBHs and the
fact that the NS radius depends weekly on their mass for hadronic equations of state (EoSs). 
A population of low-mass binaries of mostly BNSs and a few NSBHs will 
lead to inferred radii that are either approximately common (the BNSs) or outliers
(the NSBHs).
We show that BNSs and NSBHs can be identified within such a mixed population based on 
their inferred radii with high confidence, allowing us to
estimate the rate of low-mass NSBHs and achieve 
an unbiased measurement of NS radii.





§ METHOD AND RESULTS



Consider a low-mass binary with estimated component masses in the range (0.5,2)M_⊙, consistent 
with known NS masses and GW170817 <cit.>.
 In this mass range and for hadronic 
 EoSs that can support at least 2M_⊙ NSs <cit.>, 
the NS radius is expected to be constant to within a few hundred meters <cit.>. 
If the system is a BNS, then we can infer this almost-common radius, but for a misidentified
 NSBH any radius estimate will be biased. 

To quantify the bias we assume that the first binary component is a NS (the
presence of which can be confirmed by detection of an EM counterpart or tidal effects) with mass m_1 and  tidal deformability Λ_1, 
while the second component could be either a NS or a BH with mass m_2 and 
tidal deformability Λ_2 (Λ_2=0 for BHs).
In either scenario, the leading order tidal effects will be
encoded in the GW phase through


    Λ̃   ≡16/13(m_1+12m_2)m_1^4 Λ_1 +(m_2+12m_1)m_2^4 Λ_2  /(m_1+m_2)^5.


A GW analysis estimates Λ̃_est=Λ̃ if the source is a BNS, 
or Λ̃_est=Λ̃(Λ_2=0) if it is a NSBH and an error.

The NS radius is then inferred from Λ̃_est with use of two relations that do 
not sensitively depend on the EoS. 
The first relates the NS compactness to the
tidal deformability C=C(Λ), and can be used to obtain the radius from the tidal 
deformability and the mass, R=m/C(Λ) <cit.>. 
This relation holds for any NS, regardless of whether it is part of a NSBH or a BNS. 
The second relation applies to BNSs only and it relates the individual tidal 
deformabilities of the two binary components given their 
mass ratio <cit.>.

 
Working under the assumption that the binary is a BNS (a common assumption for GW170817), we
use the two EoS-insensitive relations to obtain R_ (R_), 
the radius estimate if the signal is emitted by a BNS (NSBH). The former is close to the
correct NS radius R_, while the latter is biased.
The difference between the two depends on R_ and the 
masses of the stars 


    R_-R_≡Δ R(R_,m_1,m_2)>0,


and it is plotted in Fig. 1 of <cit.>. 
The difference is smaller for larger m_2: the tidal deformability is a steeply decreasing function of the mass and almost negligible for a 2M_⊙ NS. 
Misinterpreting a heavy BH for a NS induces almost 
negligible error in the radius estimate, but NSBHs with 1M_⊙ BH result in a heavily biased radius estimate.






 §.§ Simulation of a population


Now consider a population of N low-mass binaries comprised mostly of BNSs, 
but possibly contaminated by a few 
NSBHs. Information from the BNSs will result in an unbiased estimate of the true NS radius R_ NS, 
while the corresponding radius estimate from the NSBHs will be biased by 
Δ R(R_,m_1,m_2).
To simulate such a population we assume that the signal-to-noise ratio 
(SNR), ρ, of each event follows the power-law 
distribution 3ρ_ th/ρ^4 <cit.>, 
where ρ_ th≡ 12 is the network SNR detection threshold. 
This SNR distribution is a reasonable choice since the (0.5-2)M_⊙ detectable 
binaries will be relatively local 
(redshift less than 0.1) with current GW detectors <cit.>.
 
We draw NS and BH masses from a uniform distribution in (0.5,2)M_⊙
and set all NS radii to R_=12km, consistent with the median radius 
measurement of <cit.>. 
The inferred radius for each event i has a standard deviation σ_R_i 
which is set to ∼0.75km at ρ=33, consistent with 
GW170817 <cit.> and scales inversely
with the SNR of the event. 
The likelihoood for the inferred radius of each event is then approximated with a normal distribution centered at 
R_i+N(0,σ_R_i) and with a standard deviation σ_R_i,
where R_i=R_ if the event is a BNS, or R_i=R_-Δ R(R_,m_1i,m_2i) 
if it is a NSBH. 
The additional scatter in the mean of the likelihood is caused by the random 
instance of detector noise. We approximate the likelihood for the component masses similarly, assuming 
a standard deviation of σ_m_i=0.1 M_⊙ at ρ=33 <cit.>.




 §.§ Special Event Analysis


Given the above population and corresponding radius measurements we first study whether 
we can determine the nature of individual events.
Our method is based on the fact that the inferred radii from the BNSs will be 
consistent with R_ NS, while the NSBHs result in a biased radius whose value depends on the 
component masses. 

We divide the N detections into two groups: 
a special event whose type we want to determine and the remaining N-1 detections. 
We compute the Bayes Factor (BF) that the special event is a BNS compared to a NSBH


    BF=∫ p(R'|ℋ_) L_s(d|R') dR'/∫ p(R'|ℋ_) L_s(d|R') dR',
 

where L_s(d|R') is the radius likelihood for the special event given the GW data d and 
p(R'|ℋ_) or p(R'|ℋ_) is the prior assuming the event is a BNS 
or NSBH respectively. If BF>1(<1), the GW data are more consistent with the event being a BNS (NSBH). 

The radius likelihood for the special event is computed as detailed above, while the priors are computed 
by making use of the remaining N-1 events.
We multiply the radius likelihoods for the N-1 detections and obtain the combined 
likelihood f(d|R). 
Assuming a low ratio of NSBHs to BNSs, or equivalently that the N-1 events are mostly
BNSs, f(d|R) will be consistent with R_.
Assuming a flat prior on the radius, the appropriately-normalized combined likelihood can be 
interpreted as the prior probability on
the radius for a BNS event not belonging in the N-1 detections, for example our special event: 
p(R|ℋ_)=f(d|R). If the special event is a NSBH, then the prior
can be computed again using f(d|R) and shifting it by the expected radius bias


    p(R|ℋ_)=∫   p_s(m_1,m_2|d)×
       f(d|R+Δ R(R,m_1,m_2))   dm_1 dm_2,


where p_s(m_1,m_2|d) is the posterior of the two component masses of the special event. 



We apply this method to simulated events.
We consider 1500 populations, compute the BF for each special event, and from those the probability
of correct identification. 
We find that we can correctly identify the binary type if the special event is selected wisely. 
In Fig. <ref> we consider the highest-SNR event
as this event would have small uncertainty in radius and mass. 
We find that the highest-SNR event is 
correctly identified 80% of the time after ∼ 5 events
if 20% of them are NSBHs. 
The overall probability of correct identification reaches 90% after ∼ 40 events. 
For larger ratios of NSBHs to BNSs, the NS radius prior might 
not represent the true radius. 
Such a biased measurement lowers the probability of correct identification. However, even if 
half the events are NSBHs, the probability of correctly classifying 
the highest-SNR event is ∼ 70% after about 10 detections.  


























 §.§ Hierarchical Mixture Model


The single event analysis allows for a high confidence identification/exclusion of NSBHs with a small number 
of events, however the analysis is only for identification purpose. In order to further 
measure the ratio of NSBHs to BNSs in a population and infer the NS radius we employ a hierarchical approach <cit.>.
The inferred radii follow a common underlying distribution which we model with a mixture model with two gaussian 
components and the likelihood


    L∼ (1-A) N(R_1,α_1)+ A N(R_2,α_2).


The first gaussian component models 
 the BNSs with a common radius R_1, while the second gaussian 
component models the NSBHs. We use a prior on R_1 that is uniform in [10-14] km; for R_2 we use a uniform prior in 
[R_1-10,R_1-3] since the inferred radii from NSBHs are smaller than the corresponding radii
from BNSs. The parameter A is the
ratio of NSBHs to BNSs so we use a uniform prior in [0,1]. We assume that the rate ratio 
does not evolve with redshift, a reasonable assumption for low-mass binaries detected by second generation  detectors.

The scatter α_1 in the radii of the BNSs
is caused by the detector noise realization. 
To find a suitable prior for α_1
we analyze BNS-only populations and find that the posterior for α_1 can be
approximated by a lognormal distribution with a mean of 0.8/√(N) km 
and a standard deviation of 1 km. 
The scatter in the NSBH radii α_2
is a combination of detector noise and the fact that the inferred radius from NSBHs depends on 
the bodies' masses. For lack of knowledge of the NSBH mass distribution we simply use a
wide prior for α_2: a lognormal distribution with a mean of √(3) km 
and a standard deviation of 1 km. We have verified that all prior bounds do
 not affect the resulting posteriors, with the obvious exception of A.



We simulate populations of low-mass detections that are potentially contaminated by NSBHs 
and compute the posterior of the 5 parameters of the hierarchical mixture model, Eq. (<ref>).
This method can correct the bias in the NS radius estimate even
if the population includes NSBHs as we show in Fig. <ref> which plots the posterior
for R_1 with and without (setting A=0) the mixture model for different values of the ratio of NSBHs to BNSs. 
In all cases
the mixture model is able to separate the detected events well-enough into BNSs and outliers such 
that it leads to a correct estimate of the true BNS radius. 





Besides a corrected measurement of the NS radius, we also obtain an estimate
for A, the ratio of NSBHs to BNSs. In Fig. <ref> we plot credible
intervals for A as a function of the number of events, averaged over 200 populations. We find that if no low-mass NSBHs exist we put an upper limit on their relative abundance of 3% at the 90% level with 100 detections. If, on the other hand, low-mass NSBHs do exist we can constrain their abundance to within 0.16(0.11)[0.08] at the 68% level with 100 detections if the true ratio is 0.3(0.2)[0.1].








§ DISCUSSION



We present a method to identify NSBHs in a population of low-mass events, 
measure their relative abundance, and measure the NS radius.
We find that we can correctly classify the loudest events with only a handful of detections and 
measure the ratio of NSBHs to BNSs with a few dozens of events.
In fact, the combined merger rate of GW170817 and GW190425 is 1090_-800^+1720 Gpc^-3yr^-1 <cit.>. 

This applies to any merger in this mass range, be it a 
BNS or a NSBH, suggesting a few to many tens of relevant detections
in the upcoming observing runs <cit.>. 
We therefore expect the identification of a BNS or a NSBH with GWs alone in the near future
and a measurement of their rate ratio with a few years of data. 
We emphasize that we do not use information from the post-inspiral phase of the binary, or
rely on EM counterparts to the mergers. 

Our approach treats NSBHs as outliers in a population so its performance is degraded if the
fraction of NSBHs is high.
However, we show that the probability of correct identification of the event with the largest
SNR reaches 70% after 
10 detections even if 50% of the low-mass mergers are NSBHs. 
Similarly, we find that our ratio posteriors in Fig. <ref> are systematically shifted to lower 
values of A as A increases. Despite that, we can recover the rate ratio 
at the 1σ level for a ratio up to at least 30%. Moreover, we obtain an unbiased radius estimate
event for A=30% as our approach is based on identifying radius outliers; any potentially misidentified NSBH 
will have an inferred radius consistent with the BNSs and will thus not bias the radius estimate.
	
For our simulations we assumed a true NS radius of 12 km. A stiffer EoS, a heavier BH,
or a lighter NS will lead to a larger bias in the measured NS radius <cit.> and make
 classification and measurement of the ratio A easier. 
We also assume that the NS and BH masses are distributed uniformly in (0.5,2)M_⊙. 
If the NS mass distribution instead favors heavy stars while most BHs are lighter, 
both classification and the ratio measurement will improve. We expect the contrary if
low-mass BHs have masses around 2M_⊙.
 
One caveat is that our analysis is formulated in terms of the NS radius
and the assumption that it is approximately constant for all BNSs, at least to within statistical errors.
This is reasonable for hadronic EoSs, but it is not expected to hold for EoSs with phase
transitions to quark matter <cit.>. 
We do not consider this a limitation as our analysis can also be formulated
in terms of a quantity that it truly universal for all NSs: the EoS itself. 
In fact, the radius is correlated with the pressure at twice 
the nuclear density <cit.>, suggesting that our arguments can be 
applied to the EoS directly. 
Specifically, a population of BNSs will yield an ever-improving 
measurement of the common EoS, while a misidentified NSBH will result in
 an EoS that is different than the population. 
 In addition, the tidal deformability of neutron stars with exotic matter are likely to be similar to 
 stars with hadronic matter, so we do not expect a confusion between NSBHs and BNSs with exotic matter.

Other systematic errors in the analysis might affect the inferred radius itself. 
We test this by artificially 
widening the radius likelihood by 500m, and find that 
the probability of correct identification is reduced by just ∼ 5%. 
In reality, Ref. <cit.> argued that systematic errors are small even for a loud event
like GW170817 and they are likely to remain smaller than statistical errors until we detect a 
signal with an SNR of about 100 <cit.>. 


As a final note, seeking outliers in a population can also be used to identify
exotic systems, such as binaries with at least one quark star, or hybrid binaries where one (or both)
stars have undergone a phase transition.







§ ACKNOWLEDGEMENTS


We acknowledge valuable discussions with Will Farr, Ramesh Narayan, Yu-Dai Tsai, Salvatore Vitale. 
We thank Reed Essick and Francesco Pannarale for useful comments.
H.-Y.C. was supported by by the Black Hole Initiative at Harvard University, which is funded by grants from the John Templeton Foundation and the Gordon and Betty Moore Foundation to Harvard University.

The Flatiron Institute is supported by the Simons Foundation.

Plots in this manuscript have been made with matplotlib <cit.> and we have used stan <cit.> to sample the mixture model.




