Distinguishing binary neutron star from neutron star-black hole mergers with gravitational waves

TL;DR

This paper introduces a gravitational-wave based method to distinguish neutron star-black hole mergers from binary neutron star mergers without relying on electromagnetic signals or post-inspiral data, using population analysis.

Contribution

The authors develop a novel population-based approach that identifies mixed binaries by detecting outliers in neutron star property measurements from gravitational waves.

Findings

As few as 5 low-mass events can identify a mixed binary with 80% confidence.

With 100 events, the method can constrain the fraction of mixed binaries to about 0.1.

The approach does not require electromagnetic counterparts or post-inspiral information.

Abstract

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 $\sim 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 $\sim 0.1$ at the $68\%$ credible level with 100 events.

Figures

Equations

$$\tilde{\Lambda }$$

$$R_{\text{mboxBNS}}-R_{\text{mboxNSBH}}\equiv \Delta R(R_{\text{mbox}\mathrm{NS}},m_1,m_2)>0,$$

$$\mathrm{BF}=\frac{\int p(R^{\prime }|{\mathcal{H}}_{\text{mbox}\mathrm{BNS}}){\mathcal{L}}_s(d|R^{\prime })d{R}^{\prime }}{\int p(R^{\prime }|{\mathcal{H}}_{\text{mbox}\mathrm{NSBH}}){\mathcal{L}}_s(d|R^{\prime })d{R}^{\prime }},$$

$$p(R|{\mathcal{H}}_{\text{mbox}\mathrm{NSBH}})=\int$$

$$f(d|R+\Delta R(R,m_1,m_2))d{m}_{1}d{m}_2,$$

$$\mathcal{L}\sim (1-A)\mathcal{N}(R_1,{\alpha }_1)+A\mathcal{N}(R_2,{\alpha }_2).$$

Full text

1. Introduction

The gravitational-wave (GW) event GW170817 detected by Advanced LIGO (Aasi et al., 2015) and Virgo (Acernese et al., 2015) is consistent with the merger of two neutron stars (BNS) (Abbott et al., 2017b). 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 (Abbott et al., 2018c). Arriving after the GW signal, the electromagnetic (EM) counterparts GRB 170817A (Abbott et al., 2017a; Goldstein et al., 2017; Savchenko et al., 2017) and kilonova AT 2017gfo, e.g. Coulter et al. (2017); Soares-Santos et al. (2017), 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) (Abbott et al., 2018c; Hinderer et al., 2018; Coughlin & Dietrich, 2019).

X-ray binaries suggest a lack of BHs with mass below $5M_{\odot}$ (Bailyn et al., 1998; Farr et al., 2011; Kreidberg et al., 2012), but the origin of this mass gap between BHs and NSs is not fully understood (Kreidberg et al., 2012; Belczynski et al., 2012). Recently, the discovery of a $3.3^{+2.8}_{-0.7}M_{\odot}$ unseen companion of the giant star 2MASS J05215658+4359220 further challenged the existence of the mass gap (Thompson et al., 2019). Scenarios for the production of low-mass BHs include primordial density fluctuations (Carr et al., 2016), slow supernova explosions (Belczynski et al., 2012), mergers of NSs (Faber & Rasio, 2012), and interactions of dark matter and NSs (Bramante et al., 2018). Low-mass binary mergers can potentially help study the black hole mass distribution (Abbott et al., 2018a), but probing the existence of objects in the mass gap is challenging (Littenberg et al., 2015; Mandel et al., 2015).

As already noted in Abbott et al. (2017b), though, constraining the component masses in $\sim(0.5-2)M_{\odot}$ 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 (Flanagan & Hinderer, 2008; Hinderer, 2008). 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, $\tilde{\Lambda}$ (Wade et al., 2014). For an NSBH that is not particularly loud (signal-to-noise ratio less than $\simeq 30$), the tidal deformability is generally difficult to measure  (Pannarale et al., 2011; Lackey et al., 2012, 2014; Kumar et al., 2017; Thompson et al., 2020).

GW170817 data place a lower limit on $\tilde{\Lambda}$ subject to the assumption of small spins (Abbott et al., 2018c); the data are nonetheless consistent with a highly spinning BH binary. At the same time, a nonzero $\tilde{\Lambda}$ 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 (Hinderer et al., 2018; Coughlin & Dietrich, 2019). Similar analyses for near-future detections are subject to the availability and interpretation of an EM counterpart, while post-merger information (Abbott et al., 2017c; Chatziioannou et al., 2017; Torres-Rivas et al., 2019) or evidence for disruption (Pannarale et al., 2015) 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 (Yang et al., 2018). 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 $\sim 2M_{\odot}$ BH for a NS induces a negligible error, while misidentifying a $\sim 1M_{\odot}$ 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.

2. Method and Results

Consider a low-mass binary with estimated component masses in the range $(0.5,2)M_{\odot}$, consistent with known NS masses and GW170817 (Abbott et al., 2017b, 2018c, 2018b). In this mass range and for hadronic EoSs that can support at least $2M_{\odot}$ NSs (Antoniadis et al., 2013), the NS radius is expected to be constant to within a few hundred meters (Özel & Freire, 2016). 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 $\Lambda_{1}$, while the second component could be either a NS or a BH with mass $m_{2}$ and tidal deformability $\Lambda_{2}$ ($\Lambda_{2}=0$ for BHs). In either scenario, the leading order tidal effects will be encoded in the GW phase through

[TABLE]

A GW analysis estimates $\tilde{\Lambda}_{\textrm{est}}=\tilde{\Lambda}$ if the source is a BNS, or $\tilde{\Lambda}_{\textrm{est}}=\tilde{\Lambda}(\Lambda_{2}=0)$ if it is a NSBH and an error.

The NS radius is then inferred from $\tilde{\Lambda}_{\textrm{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(\Lambda)$, and can be used to obtain the radius from the tidal deformability and the mass, $R=m/C(\Lambda)$ (Maselli et al., 2013; Yagi & Yunes, 2017). 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 (Yagi & Yunes, 2016; Chatziioannou et al., 2018).

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_{\textrm{\mbox{\tiny BNS}}}$ ($R_{\textrm{\mbox{\tiny NSBH}}}$), the radius estimate if the signal is emitted by a BNS (NSBH). The former is close to the correct NS radius $R_{\rm\mbox{\tiny NS}}$, while the latter is biased. The difference between the two depends on $R_{\rm\mbox{\tiny NS}}$ and the masses of the stars

[TABLE]

and it is plotted in Fig. 1 of Yang et al. (2018). 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_{\odot}$ NS. Misinterpreting a heavy BH for a NS induces almost negligible error in the radius estimate, but NSBHs with $1M_{\odot}$ BH result in a heavily biased radius estimate.

2.1. 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_{\rm NS}$, while the corresponding radius estimate from the NSBHs will be biased by $\Delta R(R_{\rm\mbox{\tiny NS}},m_{1},m_{2})$. To simulate such a population we assume that the signal-to-noise ratio (SNR), $\rho$, of each event follows the power-law distribution $3\rho_{\rm th}/\rho^{4}$ (Schutz, 2011; Chen & Holz, 2014), where $\rho_{\rm th}\equiv 12$ is the network SNR detection threshold. This SNR distribution is a reasonable choice since the $(0.5-2)M_{\odot}$ detectable binaries will be relatively local (redshift less than 0.1) with current GW detectors (Abbott et al., 2013).

We draw NS and BH masses from a uniform distribution in $(0.5,2)M_{\odot}$ and set all NS radii to $R_{\rm\mbox{\tiny NS}}=12$km, consistent with the median radius measurement of Abbott et al. (2018b). The inferred radius for each event $i$ has a standard deviation $\sigma_{R_{i}}$ which is set to $\sim$0.75km at $\rho=33$, consistent with GW170817 (Abbott et al., 2018b) 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}+{\cal{N}}(0,\sigma_{R_{i}})$ and with a standard deviation $\sigma_{R_{i}}$, where $R_{i}=R_{\rm\mbox{\tiny NS}}$ if the event is a BNS, or $R_{i}=R_{\rm\mbox{\tiny NS}}-\Delta R(R_{\rm\mbox{\tiny NS}},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 $\sigma_{m_{i}}=0.1M_{\odot}$ at $\rho=33$ (Abbott et al., 2018b).

2.2. 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_{\rm 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

[TABLE]

where ${\cal{L}}_{s}(d|R^{\prime})$ is the radius likelihood for the special event given the GW data $d$ and $p(R^{\prime}|\mathscr{H}_{\rm\mbox{\tiny BNS}})$ or $p(R^{\prime}|\mathscr{H}_{\rm\mbox{\tiny NSBH}})$ is the prior assuming the event is a BNS or NSBH respectively. If ${\rm 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_{\rm\mbox{\tiny NS}}$. 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|\mathscr{H}_{\rm\mbox{\tiny BNS}})=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

[TABLE]

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. 1 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 $\sim 5$ events if $20\%$ of them are NSBHs. The overall probability of correct identification reaches $90\%$ after $\sim 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 $\sim 70\%$ after about 10 detections.

2.3. 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 (Loredo, 2004) . The inferred radii follow a common underlying distribution which we model with a mixture model with two gaussian components and the likelihood

[TABLE]

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 $\alpha_{1}$ in the radii of the BNSs is caused by the detector noise realization. To find a suitable prior for $\alpha_{1}$ we analyze BNS-only populations and find that the posterior for $\alpha_{1}$ can be approximated by a lognormal distribution with a mean of $0.8/\sqrt{N}$ km and a standard deviation of $1$ km. The scatter in the NSBH radii $\alpha_{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 $\alpha_{2}$: a lognormal distribution with a mean of $\sqrt{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. (5). This method can correct the bias in the NS radius estimate even if the population includes NSBHs as we show in Fig. 2 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. 3 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].

3. 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* (Abbott et al., 2020). 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 (Abbott et al., 2013). 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. 3 are systematically shifted to lower values of $A$ as $A$ increases. Despite that, we can recover the rate ratio at the 1$\sigma$ 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 (Yang et al., 2018) 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_{\odot}$. 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_{\odot}$.

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 (Han & Steiner, 2018). 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 (Özel & Freire, 2016), 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 $\sim 5\%$. In reality, Ref. Abbott et al. (2018c) 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 (Dudi et al., 2018).

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.

4. 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 (Hunter, 2007) and we have used stan (Carpenter et al., 2017) to sample the mixture model.