Using gravitational waves to distinguish between neutron stars and black holes in compact binary mergers

In August 2017, the first detection of a binary neutron star merger, GW170817, made it possible to study neutron stars in compact binary systems using gravitational waves. Despite being the loudest gravitational wave event detected to date (in terms of signal-to-noise ratio), it was not possible to unequivocally determine that GW170817 was caused by the merger of two neutron stars instead of two black holes from the gravitational-wave data alone. That distinction was primarily due to the accompanying electromagnetic counterpart. This raises the question: under what circumstances can gravitational-wave data alone, in the absence of an electromagnetic signal, be used to distinguish between different types of mergers? Here, we study whether a neutron star-black hole binary merger can be distinguished from a binary black hole merger using gravitational-wave data alone. We build on earlier results using chiral effective field theory to explore whether the data from LIGO and Virgo, LIGO A+, LIGO Voyager, the Einstein Telescope, or Cosmic Explorer could lead to such a distinction. The results suggest that the present LIGO-Virgo detector network will most likely be unable to distinguish between these systems even with the planned near-term upgrades. However, given an event with favorable parameters, third-generation instruments such as Cosmic Explorer will be capable of making this distinction. This result further strengthens the science case for third-generation detectors.


INTRODUCTION
Neutron stars are unique laboratories for studying ultra-dense, relativistic matter. Multimessenger observations of neutron star mergers provide unique opportunities to extract relevant physical information (such as compactness) from them. Measurements of neutron star compactness and radii are vital to constraining the equation of state of ultra-dense matter Lattimer & Prakash (2001). In addition to gravitational wave observations of the neutron star mergers GW170817 and GW190425 Abbott et al. (2017, 2020, X-ray observations of accreting neutron starsÖzel & Freire (2016); Watts et al. (2016) have placed constraints on neutron star mass and radii. Of these, the recent results from NICER are especially promising Bogdanov et al. (2019a,b); Raaijmakers et al. (2020).
The observation of GW170817 and its electromagnetic counterpart led to several important advances. The detection of the electromagnetic counterpart was possible because the LIGO-Virgo observation constrained the sky location of the event to 28 deg 2 . It was the detection of gamma-ray burst GRB170817A 1.7 seconds after GW170817 that provided the initial evidence that this event contained neutron star matter. Transient electromagnetic follow-ups Soares-Santos et al.  Abbott et al. (2018). Later work constrained the properties of this event further M 1 = 1.45 +0.08 −0.06 , R 1 = 12.36 +0.52 −0.38 and M 2 = 1.28 +0.05 −0.06 , R 2 = 12.32 +0.66 −0.43 Fasano et al. (2019). Combining gravitational wave observations GW170817 and GW190425 with NICER results led to constraints on the radius of a 1.4M of 12.33 +0.76 −0.81 km and 12.18 +0.56 −0.79 km Raaijmakers et al. (2021). Though GW170817 led to new constraints on the radii and tidal deformabilities of neutron stars, it alone was not sufficient to determine that the event was a binary neutron star rather than a binary black hole. The evidence that this was a binary neutron star merger came from observations of the electromagnetic counterpart. In future observations, we will likely not be in the fortuitous position of having a clear electromagnetic counterpart. To date, LIGO-Virgo has detected two neutron star-black hole binaries, neither of which had an electromagnetic counterpart Abbott et al. (2021). Furthermore, if the mass of any of the binary components happen to lie within the mass gap, gravitational waves are the most promising avenue by which to determine whether the object is a black hole or a neutron star. This leads to the questions: under what conditions can a gravitational wave signal alone differentiate between a binary neutron star and a binary black hole? Can a neutron star-black hole binary be differentiated from a binary black hole by gravitational wave observations alone? This work addresses the second of these questions for current and future gravitational wave detectors. Current detectors may not be able to successfully differentiate between neutron star-black hole binaries and binary black holes, making future detectors vitally important. The importance of future detectors for studying neutron stars in binary neutron star mergers was shown in a recent paper Pacilio et al. (2021).
In addition to the current LIGO-Virgo detectors, we consider LIGO A+, LIGO Voyager, the Einstein Telescope, and Cosmic Explorer. The plans for LIGO A+ aim to improve the detection range of binary neutron stars at 1.4 M by a factor of 1.9 Barsotti et al. (2018). These improvements to LIGO may occur as soon as three years from now. Further plans exist for LIGO Voyager, which will further increase detector sensitivity McClelland et al. (2016). Power spectral density curves for the design sensitivity of these two detectors are publicly available Evans et al. (2020) and are used in our analysis. Beyond LIGO A+ and LIGO Voyager, there are plans for third-generation (3G) detectors such as the Einstein Telescope (ET) and Cosmic Explorer. We consider the Einstein Telescope and both Cosmic Explorer's first run (CE1) expected to take place in the 2030s and its second run (CE2) which is planned for the 2040s. Cosmic Explorer is expected to vastly increase the number of neutron stars detected by expanding the redshift horizon for binary neutron star detections out to 3.1 in the first run. With predicted signal-to-noise ratios going up by an order of magnitude for nearby sources, third-generation detectors will significantly improve our tidal deformability measurements Reitze et al. (2019).
We use standard Bayesian model selection tools in our analysis. The evidences are calculated using the dynamic nested sampling package DYNESTY Speagle (2020) accessed through the PyCBC toolkit Biwer et al. (2019). In this analysis, we employ neutron star equations of state derived from chiral effective field theory, a theory that uses an effective description of nuclear matter in terms of nucleons and pions Weinberg (1990Weinberg ( , 1991; Machleidt & Entem (2011);Epelbaum et al. (2009). The chiral effective field theory framework not only leads to equations of state consistent with all symmetries of the strong interactions and known experimental constraints, but it also provides reliable uncertainty estimates. We use the same subset of the equations of state employed successfully in Capano et al. (2020) to improve constraints on neutron star radii.
We show that, at least for the proposed gravitational wave detectors within the next decade (namely LIGO A+ and Voyager), it is very unlikely that gravitational wave observations alone will be able to distinguish neutron star-black hole binaries from binary black holes. Third-generation gravitational wave detectors will be required for this purpose. Sec. 2 details our model selection procedure, Sec. 3 presents the main results, and Sec. 4 concludes with a discussion of the implication of these results.

METHODS
Consider a network of gravitational wave detectors, and let d i (t) denote the gravitational wave strain time series data in the i th detector as a function of time t. The collection of all time series data in the network will be denoted d. The data is the sum of detector noise n i (t) and a possible astrophysical signal h(t), which depends on certain parameters which we collectively denote ϑ: The central goal of a Bayesian analysis is to calculate the posterior probability distributions p( ϑ| d(t)) of the parameters ϑ. The basis of this is Bayes' Theorem The fourteen parameters appearing in ϑ are discussed below. The prior, p( ϑ|H), represents the knowledge that we have about the parameters before considering the data. The likelihood function, p( d(t)| ϑ, H), is the probability of obtaining the observation d(t) given a waveform H with parameters ϑ.
In order to obtain a posterior distribution on one or a few parameters, we marginalize over the other parameters by integrating p( d(t)| ϑ, H)p( ϑ|H). Marginalizing over all parameters yields the evidence. Comparing the evidence (p( d(t)|H)) of two different models (H A and H B ) gives the Bayes factor, This number indicates how much the data supports one model over the other. When B > 1, H A is favored over H B ; the larger B is, the more H A is favored. In this study, the Bayes factors express how much A, the neutron star-black hole model, is favored over B, the binary black hole model. We measure evidences using the dynamic nested sampling package DYNESTY Speagle (2020). To crosscheck our results, we analyze a subset of our signals using a parallel-tempered version of the emcee Markov Chain Monte Carlo sampler Vousden et al. (2015); Foreman-Mackey et al. (2013). The resulting posteriors were consistent with those generated by DYNESTY. We generate simulated gravitational waves from neutron star-black hole binary (NSBH) mergers and add these to simulated Gaussian noise colored by the target detector configuration's power spectral density (PSD). Gravitational waves from neutron star-black hole mergers depend on multiple variables ϑ. The most relevant parameters for this work are the component masses m 1,2 and the tidal deformabilities Λ 1,2 , defined as Here R 1,2 are the radii of the individual stars, and k 2 is the tidal Love number, which is determined from the equation of state and the mass. The leading order effect of Λ 1,2 on the waveform is through the combined tidal deformability parameter where we define the mass ratio q = m 2 /m 1 ≥ 1. The tidal deformability is the primary means to distinguish black holes from neutron stars using gravitational waves and infer the equation of state of neutron stars. By definition, a black hole has zero tidal deformability, while larger values of Λ correspond to stiffer equations of state. As the binary inspirals, merges, and then settles into a stable black hole, it emits gravitational waves. Gravitational waves have two polarizations, denoted as h +,× . The intrinsic parameters affect the phase evolution of the gravitational waves. Some parameters, such as the chirp mass and symmetric mass ratio, also affect the amplitude of the gravitational wave. The symmetric mass ratio ν and chirp mass M are defined respectively as In the source frame, say one aligned with the source axis, h +,× depend on the direction to the detector; equivalently, in a geocentric frame, h +,× depend on the orientation of the source. Furthermore, the detectors do not detect h + and h × directly, they detect the gravitational wave strain: where F + and F × are functions of the angles defining the location of the source. These angles are typically expressed as sky location (right ascension α, declination δ) and polarization angle ψ. Additionally, the amplitude depends on the inclination angle ι and the luminosity distance of the source. These extrinsic variables affect only the amplitude of the gravitational waveform. The last variable that defines the detected strain is the detection time t c (which determines the detector position and orientation). We employ the PyCBC toolkit Biwer et al. (2019) to generate the gravitational waveforms. This requires a specification of the parameters ϑ and a waveform approximant. In the data analysis, colored Gaussian noise is added to the generated waveform using the detector power spectral density curve. We use a waveform approximant that combines inspiral, merger, and ringdown portions of the signal and has been calibrated to numerical relativity results (see e.g. Ajith et al. (2008) Thompson et al. (2020) were considered as well; however, we found SEOB-NRv4 ROM NRTidalv2 NSBH to be the best choice for this analysis. We note here that neutron star-black hole systems present a considerable challenge for existing signal models and significant uncertainties remain. This is especially true at the high signal-to-noise ratios possible for third-generation detectors. For this reason, we compare the results for all three waveforms.
We set the neutron star mass to the standard 1.4M and vary the mass of the black hole between (5, 10, 15, 20)M and the distance between (40,80) Mpc. For the neutron star, we choose two of the equations of state based on chiral effective field theory that were favored by parameter estimation in a previous work Capano et al. (2020). The first of these equations is the maximum For the injection, the polarization, inclination, and coalescence time are set to Ψ = π , tc = 1187008882.4434 , ι = 0.35 .
The choice of sky location and inclination is arbitrary, and the effects of choice of sky location are discussed in Section 4. For the analysis with DYNESTY, we set up the parameter estimation to be as similar to the analysis for GW170817 as possible. As was done for GW170817 Capano et al. (2020); De et al. (2018), we fix the sky location and distance. While it is unlikely that the sky location of a detected neutron star-black hole system will be known to such accuracy, fixing the sky location in the analysis significantly reduces computation time and does not effect the resulting Bayes factors. To confirm this, we performed a series of parameter estimation runs where the sky location was a variable parameter and found the Bayes factors to be completely consistent. The variable parameters in our parameter estimation were the individual masses, spins, coalescence time, inclination, and polarization. The prior for the neutron star mass object was uniform on [1M ,2M ] and for the black hole it was uniform on  2017)). We constrained the inclination and polarization angles to be between 0 and 2π rad, and the coalescence time was assumed to be in the range tc ± 0.1 sec.
The tidal deformability parameter estimation is what differs between our two models. To test the binary black hole hypothesis, the tidal deformability of both objects is set to 0 in the parameter estimation. We looked at two cases for the neutron star-black hole model parameter estimation. In one case, we sampled over the equation of state for the neutron star mass object. The equation of state has a uniform prior in radius at 1.4M , and there are 2,000 equations in the prior. The equation selected by the sampler was then used to calculate the tidal deformability given m 1 and m 2 . In the other case, we set the equation of state as a static variable in the parameter estimation. The reason for this is that, while the nuclear equation of state is currently not well constrained, it is expected that experiments such as NICER will significantly improve our knowledge over the next decade. To take this into account, we consider the extreme case: the one in which the equation of state is known exactly and is thus fixed in the parameter estimation.

RESULTS
To determine if gravitational waves can distinguish between neutron star-black hole binaries and binary black holes, we look at the natural log of the Bayes factor (ln B) between two models. There is much debate on what constitutes evidence, strong evidence, decisive evidence, and so on. Commonly cited statistics papers such as Kass & Raferty (1995) state that log 10 B ≥ 2 (ln B ≥ 5) can be considered decisive evidence in favor of a model. However, this is questionable for gravitational wave model selection because of the high dimensionality, complexity and several degeneracies of the parameter space (which are not yet fully understood). Additionally, using different sampler settings and different noise realizations can lead to variations in ln B of about ±2 at the 1σ level when ln B ≈ 5, and around ±4 when ln B ≈ 10. Taking account these uncertainties, we have decided to require a higher threshold thereby ensuring that our conclusions remain conservative regarding the capabilities of the gravitational wave detectors that we consider. We require ln B ≥ 10 for decisive evidence. The errors quoted in this paper are based on the standard deviation of ln B across instances of the same injection parameters but with different noise realizations. Except for the specific case of m BH = 5M , the errors for the current LIGO-Virgo detector network, LIGO A+, and LIGO Voyager are < 1. For the Einstein Telescope and Cosmic Explorer 1, the errors are < 2.5, and for Cosmic Explorer 2, the errors are < 4. The errors in the case of m BH = 5M are larger (for details see Table 1). The maximum error for the current LIGO-Virgo detector network is ≈ 1.0. This increases to 2.5 for LIGO A+, 4.7 for LIGO Voyager, ≈ 16 for Einstein Telescope, ≈ 19 for Cosmic Explorer 1, and ≈ 29 for Cosmic Explorer 2. The relative error decreases by nearly an order of magnitude as signal-to-noise ratio increases, i.e. from 1 for aLIGO and Virgo to 0.1 for CE2.
As mentioned earlier, we shall present results for a 1.4M neutron star with a black hole companion of mass (5, 10, 15, 20)M , and we shall take the distance to be 40 or 80 Mpc. The neutron star equation of state shall be either of the ones shown in Fig. 1. We shall consider the following detector networks: • The current LIGO-Virgo Network at design sensitivity in the zero-detuned high power configuration LIGO Scientific Collaboration ( • The second observational run of the proposed 40 km Cosmic Explorer detector, again in the "compact binary" configuration Srivastava et al. (2020).
The results of our analysis for the various combinations of masses, distances and detector network are shown in Figs. 2, 3, and 4, and in Tables 2-13. The Figs. 2 and 3 show the Bayes factors for the 5M and 10M black hole cases respectively. As expected, the 5M case leads to larger Bayes factors since the tidal effects on the neutron star are more significant. Nevertheless, for both cases, the important observation for our purposes is that the Bayes factors exceed our chosen threshold of Eq. 10 for the third-generation detectors (The Einstein Telescope and Cosmic Explorer). The Voyager results in Fig. 2 for the stiff equation of state surpass the threshold slightly. However, the variation is seen to be large. Furthermore, this is not the case for the maximum likelihood equation of state or for a 10 solar mass black hole companion. In a fine-tuned case, LIGO Voyager might be able to do this measurement. However, as it requires the event to be closer than any binary detected to date, to have a black hole companion that is smaller than any black hole observed by LIGO thus far, and to have a nuclear equation of state that is rather optimistically stiff, it is unlikely. This same conclusion is evident in Fig. 4, which shows all the combinations that we have considered: ln B > 10 almost exclusively for the third-generation detectors.
The precise numerical values for the Bayes factors are found in Tables 2-13. Looking at Tables 2 and 3, we see that for the current LIGO-Virgo detector network | ln B| < The results using the third-generation detectors (the Einstein Telescope, Cosmic Explorer 1, and Cosmic Explorer 2) are more optimistic. We finally see multiple values above the threshold of 10, though in two cases, the 1σ error falls below the cutoff. Looking at Tables 8 and 9, we see that for the Einstein Telescope, there are now three instances for both the variable and constant equation of state cases that exceed our Bayes factor threshold. In all cases, this occurs for a black hole companion of mass 5M . For the stiff equation of state, we have ln B of (81.8, 82.9) at 40Mpc and (17.  Finally, we note that current waveform models for neutron star-black hole binaries have limitations in high signal-to-noise ratio and high mass ratio regimes such as the ones explored in this paper. Developing more accurate waveform models is important for analyzing real data. However, in our studies, we inject simulated signals in noise and recover them with the same signal model. Thus it is most important that the signal model capture the same qualitative features as the true signal. We considered two neutron star-black hole models, SEOBNRv4 ROM NRTidalv2 NSBH and IMRPhenomNSBH, as well as the older IMRPhenomD NRTidal model. Figure 5 compares the results from all three waveforms for the m BH = 5M case. For LIGO-Virgo and its upgrades, all three waveforms agree at the 2σ level. For the 3G detectors, however, IMRPhenomD NRTidal gives significantly lower results than the two neutron star-black hole waveforms. We see that for all three waveforms, the results qualitatively agree: The Bayes factors for Einstein Telescope and the first and second runs of Cosmic Explorer remain comfortably above the threshold. However, a closer look reveals that IMRPhenomNSBH is unsuitable for use with 3G detectors when m BH > 10M . As the mass of the black hole companion increases, tidal effects decrease, and the gravitational waves emitted grow more sim-  ilar to those of a binary black hole system. This means that the Bayes factor of the neutron star-black hole case over the binary black hole case should approach one for large black hole masses. This is indeed the behavior observed with SEOB-NRv4 ROM NRTidalv2 NSBH and IMRPhenomD NRTidal. However, for IMRPhenomNSBH, we find that log B increases as the black hole companion mass increases from 10 to 20 M (see Figure 6). This is clearly unphysical behavior and deserves further explanation 1 . Digging still deeper, the problem turns out to be the gravitational wave amplitude; IMRPhenomNSBH uses an older ansatz for the amplitude Santamaria et al. (2010) (which was not originally intended for neutron star -black hole systems). Figure 7 compares the gravitational wave amplitude as a function of frequency for two different values of Λ, for the m BH = 20M case, for both approximants (along with the amplitude spectral density for the CE1 detector). We clearly see that while the SEOBNRv4 ROM NRTidalv2 NSBH model shows no dependence of the amplitude on Λ, which is what we expect for these high mass configurations, the IMRPhenomNSBH model shows a large dependence on Λ for the post-merger signal which is clearly unphysical. This incorrect behavior explains the effect shown in Figure 6. Due to this non-physical behavior, we choose to use SEOBNRv4 ROM NRTidalv2 NSBH for our analysis.

DISCUSSION
The results demonstrate that the current LIGO and Virgo detectors are not sufficient to differentiate between neutron star-black hole and binary black hole systems. In fact, the success of A+ and Voyager for this purpose is dubious. There were no cases for either LIGO-Virgo or LIGO A+ where the ln B exceeded our threshold, and only a fine-tuned case for LIGO Voyager. However, it is important to note that the cases with the highest ln B always occur with the stiff equation of state, the 5M black hole companion, and at 40Mpc. This is not surprising. The stiff equation of state was selected specifically for this property, and the signal-to-noise ratio at 40Mpc is higher than at 80Mpc.
It can be seen that the ability to differentiate between a neutron star-black hole system and a binary black hole system does not directly correspond to the signal-to-noise ratio. The highest ln B occurs for m BH = 5M , even though systems with m BH = 10, 15, 20M have higher signal-to-noise ratios. The tidal effects decrease as mass increases, and this effect is clearly of greater importance than the increase in signal strength. Detection of a neutron star-black hole system with a low mass ratio will almost certainly be required to give evidence of neutron star matter in the gravitational wave signal. Additionally, the nuclear equation of state itself is an important factor in how soon we will and how likely we are to distinguish a neutron star-black hole system from a binary black hole system. Finally, in this analysis, we have made a particular choice of sky location and inclination angle of the source. This orientation, corresponding to GW170817, is a favorable one. As expected, repeating the simulations with randomly chosen sky-positions generally leads to smaller Bayes factors. However, even in this case, the Bayes factors for the Einstein Telescope, Cosmic Explorer 1 and Cosmic Explorer 2 remain comfortably above the threshold, while for Voyager, the results get closer to the threshold. Our basic conclusions therefore remain unchanged.
When looking at LIGO Voyager, we saw that in one case the results were close to our cutoff. Keep in mind, however, that this analysis was done with design sensitivity curves and this result occurs only for the fine-tuned case of a very close binary with a very small black hole that has a rather stiff equation of state. Despite LIGO-Virgo's recent detection of an object in the mass gap, 5M is still on the low end of what we expect for black hole masses. Looking at Figure 3, it's evident that, when the companion mass increases to even 10M , the Bayes factor drops rapidly regardless of distance or equation of state for all detectors. If the equation of state is as soft as the analysis of GW170817 suggests, then LIGO Voyager will certainly be unable to distinguish neutron star-black hole systems from binary black holes regardless of how close or loud the signal is.
3G detectors will likely be required to obtain decisive evidence of neutron star-black hole system from gravitational wave data. We see here that the proposed designs for the Einstein Telescope and Cosmic Explorer may very well allow for these detections. Looking at Figure 2 and Table 10, we see that regardless of the nuclear equation of state, there are systems which have ln B > 10. Thus, if current analyses of GW170817 are accurate, we will be waiting until the Einstein Telescope or Cosmic Explorer for gravitational wave evidence of neutron star-black hole systems. Additionally, 3G detectors seem to be able to do this measurement at distances out 80 Mpc (stiff equation of state), which greatly expands the number of candidate systems. Even with very sensitive future detectors, the ability to distinguish neutron star-black hole systems from binary black holes is very dependent on the mass of the black hole in the binary.