Model comparison from LIGO–Virgo data on GW170817’s binary components and consequences for the merger remnant

GW170817 is the very first observation of gravitational waves originating from the coalescence of two compact objects in the mass range of neutron stars, accompanied by electromagnetic counterparts, and offers an opportunity to directly probe the internal structure of neutron stars. We perform Bayesian model selection on a wide range of theoretical predictions for the neutron star equation of state. For the binary neutron star hypothesis, we find that we cannot rule out the majority of theoretical models considered. In addition, the gravitational-wave data alone does not rule out the possibility that one or both objects were low-mass black holes. We discuss the possible outcomes in the case of a binary neutron star merger, finding that all scenarios from prompt collapse to long-lived or even stable remnants are possible. For long-lived remnants, we place an upper limit of 1.9 kHz on the rotation rate. If a black hole was formed any time after merger and the coalescing stars were slowly rotating, then the maximum baryonic mass of non-rotating neutron stars is at most 3.05 M (cid:31) , and three equations of state considered here can be ruled out. We obtain a tighter limit of 2.67 M (cid:31) for the case that the merger results in a hypermassive neutron star.

GW170817 is the very first observation of gravitational waves originating from the coalescence of two compact objects in the mass range of neutron stars, accompanied by electromagnetic counterparts, and offers an opportunity to directly probe the internal structure of neutron stars. We perform Bayesian model selection on a wide range of theoretical predictions for the neutron star equation of state. For the binary neutron star hypothesis, we find that we cannot rule out the majority of theoretical models considered. In addition, the gravitational-wave data alone does not rule out the possibility that one or both objects were low-mass black holes. We discuss the possible outcomes in the case of a binary neutron star merger, finding that all scenarios from prompt collapse to long-lived or even stable remnants are possible. For long-lived remnants, we place an upper limit of 1.9 kHz on the rotation rate. If a black hole was formed any time after merger and the coalescing stars were slowly rotating, then the maximum baryonic mass of non-rotating neutron stars is at most 3.05 M , and three equations of state considered here can be ruled out. We obtain a tighter limit of 2.67 M for the case that the merger results in a hypermassive neutron star.
Keywords: neutron stars, neutron star equation of state, gravitational wave astronomy, compact object mergers (Some figures may appear in colour only in the online journal)

Introduction
On August 17, 2017, the Advanced LIGO [1] and Advanced Virgo [2] observatories detected a gravitational-wave signal, GW170817, from a low-mass compact binary system [3]. The signal was coincident with the gamma-ray burst GRB 170817A [4] and was followed by other observations spanning the electromagnetic spectrum [5]. These electromagnetic observations, coupled with measurements of the component masses from the gravitational-wave observation [6,7], provide compelling evidence that GW170817 was produced from the merger of two neutron stars.
One of the key goals of observing such neutron star mergers is to extract information regarding the internal structure of neutron stars and to improve our understanding of their highly uncertain equation of state [8,9]. Information regarding the equation of state is encoded in the variation of the orbital evolution induced by tidal interactions and the deformation of the neutron stars under their own angular momenta [10][11][12][13][14]. Such interactions are small at low frequencies but rapidly increase in strength in the final tens of gravitational-wave cycles before merger. Additional information about the equation of state is carried in the post-merger signal, as the remnant body can be a black hole, a stable neutron star or a meta-stable neutron star that collapses after some delay [15][16][17][18][19][20][21]. However, gravitational waves from the postmerger remnant are expected at higher frequencies where the detectors are less sensitive. No such post-merger signal was observed for GW170817 [6,20,22,23].
The difficulty of computing the equation of state from first principles has led to a wide range of plausible theoretical models based on nuclear theoretical calculations. Astronomical observations of pulsars and x-ray binaries (e.g. [8,24,25]) have already been able to rule out a subset of proposed models, most notably by measurements of large neutron star masses [26]. Several studies have derived constraints on the equation of state through direct measurements of the neutron star mass and tidal deformability as inferred from GW170817 [27][28][29][30][31][32][33][34][35][36][37][38][39][40][41].
In this paper, we address whether we can rule out specific equations of state from the gravitational-wave data alone and whether we can rule out the possibility that one or both of the coalescing objects was a low-mass black hole. We opt to work with a subset of equations of state representing nuclear matter properties that vary over a wide range of possible values. The correlations between the properties of the neutron star with the underlying equation of state manifests itself in a wide range of macroscopic physical observables, such as the mass, radius or tidal deformability. In practice, we do not expect the true equation of state to be contained in any of the models used in this analysis. By comparing the Bayesian odds ratio for different models, we expect that equations of state that are most similar to the real equation of state will lead to the highest odds ratios whereas equations of state that differ significantly from the data will be down-ranked in the analysis, giving insight into the physical properties of the true equation of state.
We apply Bayesian techniques [42,43] to compare a large set of equations of state against an alternative hypothesis. The first fully Bayesian study of the feasibility of probing the neutron star equation of state was presented in [44] and expanded upon in [45]. We consider two possible baseline hypotheses: the first hypothesis is that both component objects were lowmass black holes; the second hypothesis is that GW170817 was composed of two neutron stars with an equation of state corresponding to one of the best fit equations of state in our analysis. In order to make quantitative statements about the relative likelihoods of different equations of state, we do not allow the tidal parameters to vary independently and instead enforce that each neutron star obeys a specific equation of state, similar to [38,39,[44][45][46]. This is in direct contrast to the analysis presented in [6], where the tidal parameters were allowed to vary independently implicitly allowing each neutron star to have a different equation of state. By performing parameter estimation using a fixed equation of state, we are able to conveniently incorporate a number of physical prior constraints, as was discussed in detail in [39].
We find that the equations of state yielding excessively large tidal deformabilities are disfavored by the data. For example, the MS1_PP model [47] is disfavored with a Bayes factor that is less than ∼0.005 compared to the most favored equation of state models. However, we find that we cannot comprehensively rule out the majority of equations of state used in this analysis. Instead, we find that many of the equations of state have comparable evidences that are within an order of magnitude. This includes the models with the lowest tidal deformabilities as well as the binary black hole case.
In a previous analysis of GW170817 [39], we computed radii and tidal deformabilities for GW170817 using two methods. The first method links the tidal deformabilities of the two neutron stars by assuming that they obey the same equation of state [48,49] and the second method utilized a parametrized equation of state [50,51]. The neutron star radii were found to be no larger than ∼13 km, disfavoring large radii and stiffer equations of state. In this paper, we perform a direct model comparison between equation of state models, finding results that are consistent with [39]. The equation of state models that are most preferred by the data in this analysis are the ones that predict radii close to those computed in [39].
In addition, we also explore the inferences that can be made about the fate of the remnant object for each proposed equation of state. The range of maximum neutron star masses for the equation of state models allows remnants which are indefinitely stable neutron stars, longlived supramassive neutron stars, or short-lived hypermassive neutron stars. Prompt collapse remains a possibility, though it can be ruled out for many models. In other cases, the required mass threshold from numerical relativity is not available or is insufficiently accurate to make a clear statement. For cases resulting in a supramassive neutron star, we use the mass range inferred for each of our equation of state models to derive the maximum rotation rate of the neutron star after it settles into a uniformly rotating state, obtaining an upper limit of 1.9 kHz.
Finally, we investigate constraints for the maximum mass of non-rotating neutron stars that follow from different assumptions about the type of the remnant. This is motivated by one of the possible explanations for the short gamma-ray burst models, which requires the formation of a black hole surrounded by a disk [52]. The other model assumes a rapidly rotating magnetar instead [53]. Which one is the correct description for the case at hand is an open question, but if we do assume the first one, then the presence of a black hole implies that neutron stars with the total baryonic mass inferred for the initial binary cannot be stable. The resulting limit is not very restrictive, but would rule out the three most extreme of our equations of state when assuming that the coalescing stars are slowly rotating. The bound can be further tightened by limiting the lifetime of the remnant to the observed delay of the short gamma-ray burst. The steep gradient of lifetime as function of mass near the maximum rotating neutron star mass then places the remnant mass in (or at least close to) the hypermassive range. Our limits are compatible with similar calculations carried out in [27,29,54], but are obtained using the fully consistent probability distributions computed for each equation of state model.

Gravitational-wave Bayesian model comparison
In this work our goal is to perform a direct comparison of theoretical models for the neutron star equation of state against the gravitational-wave data for GW170817. In order to do so, we employ Bayesian model selection to test a specific equation of state against a competing hypothesis. We adopt both the low-mass binary black hole hypothesis and a best fit equation of state as baseline hypotheses and assess the impact this has on the inferences that can be made from the data. Detailed introductions to Bayesian methods in gravitational wave astronomy can be found in [42,[55][56][57].
Given the observed gravitational-wave data s and any prior information I, the posterior odds in favor of a given model hypothesis M i can be expressed using the odds ratio where Z = P(s|M i , I) is the marginalized likelihood or evidence and P(M i |I) is the prior probability on the model hypothesis. If we assume that the two models are a priori equally likely, P(M 1 |I) = P(M 2 |I), then the odds ratio O 12 reduces to the Bayes factor B 12 . The evidence is calculated by marginalizing over all parameters weighted by their prior probabilities P(ξ|M i , I) where the likelihood of obtaining the data realization s given a gravitational-wave signal with physical parameters ξ is Together with a given set of prior probabilities, we have all the information needed to evaluate the evidence in equation (2). In the following subsections, we discuss the waveform models used to describe the gravitational-wave signal h(t) and their dependence on the underlying physical parameters ξ.

Prior probabilities
An important consideration when computing the Bayes factors is the choice of prior probabilities. This choice will enter the analysis in two distinct ways. First, we must make a choice for the relative prior probabilities between any two models; P(M i |I). A wide range of equations of state are still consistent with observations [58], though additional theoretical considerations and observations could allow us to restrict the range of plausible equations of state [3,6,59].
Here, neglecting such additional information, we assume that the prior probability for all equations of state is the same. As detailed above, the odds ratio then reduces to the Bayes factor between any two models. The second choice of priors that we must make is the choice of prior probabilities on the parameters that enter a specific model, P(ξ|M i , I). As the choice of priors will affect the resulting Bayes factors, we explicitly detail the choices made in our analysis below.
The choice of priors that we use in this work is broadly similar to our previous analysis of GW170817 [3,6]. We assume a priori that the source of GW170817 resides in NGC 4993, allowing us to adopt a fixed sky location. As in [6], although the distance to NGC 4993 is known, we infer the luminosity distance purely from the gravitational-wave data. The distance prior assumes that sources are uniformly distributed in volume up to a maximum distance of 75 Mpc. The coalescence time t c and phase φ c are taken to be uniformly distributed and we assume that the binary orientation and polarization angle are isotropic.
The remaining parameters that we sample are the masses and spins of the two compact objects. The choice of prior ranges for these parameters is not so easily made. Observations of neutron stars in our galaxy can give us indications of the mass and spin distributions of binary neutron star systems [8,[60][61][62][63][64][65][66], but there is some uncertainty in these measurements. Additionally, these observations would not be appropriate to use if we want to assess the possibility that the system might be two black holes.
For these reasons we consider two sets of priors in this work. The first narrow prior is motivated by observed binary neutron star systems in our galaxy [26,67]. For this prior, we assume that the prior on the component masses is Gaussian with a mean of 1.33M and a standard deviation of 0.09M [8]. The prior on the dimensionless spin magnitudes are taken to be uniformly distributed between 0 and 0.05 and the spin directions are assumed to be isotropic. This spin prior is consistent with the observed population of binary neutron stars that will merge within a Hubble time [3,68]. As many of the waveform models we consider do not allow for spin components perpendicular to the orbital angular momentum, we choose to sample the spin direction and magnitude as described above and then ignore any spin component perpendicular to the orbital angular momentum. In practice, this amounts to choosing a non-uniform distribution on the spin magnitude parallel to the orbital plane, with a peak at 0.
The second broad prior allows for a wider range of masses and spins. Here we assume that component masses are uniformly distributed between 0.7M and 3.0M with the constraint that m 2 m 1 . Priors on the chirp mass M and mass ratio q are determined from the concomitant Jacobian transformation. For computational efficiency, we impose a further cut on the chirp mass of the binary, such that we only consider points with chirp mass between 1.190 and 1.210 solar masses. However, this is wide enough to cover the full range of values supported by the data [3,6,38,39]. In the broad prior, the spin magnitudes are taken to be uniformly distributed between 0 and 0.7. This is narrower than the spin prior used in [6] which extended to 0.89. The broad prior is intended to be a conservative spin prior compatible with the fastestspinning known neutron star, which has a dimensionless spin 0.4 [69]. The choice of the limit 0.7 should not influence the results, as we obtain one-sided 95% upper limits for the dimensionless spin posterior which are much lower (below 0.35 for all models).
Non-rotating neutron stars cannot exceed a maximum mass which depends on the equation of state. If the sampler chooses a point in the parameter space where a given equation of state cannot support a non-rotating neutron star, we assume that the body is a black hole. Although rotation can increase the maximum mass, we do not allow for a scenario in which any neutron star is so massive that it requires rotation to prevent collapse. This choice is made to ensure that the prior volume is equal for all the equations of state. Our exact model assumption for a given equation of state is thus that objects below the maximum mass of non-rotating neutron stars are in fact neutron stars described by the given equation of state, otherwise light black holes. For the narrow prior, the distinction is irrelevant because we find no posterior support above the maximum mass for any equation of state. For the broad prior on the other hand, the posterior has significant support for mixed binaries for some equation of state models, as is further discussed in section 3.
There is also an upper bound on the dimensionless spin possible for uniformly rotating neutron stars, which depends on the equation of state. The maximum spin covered by the broad prior is below the maximum possible spin for some models considered here (compare section 2.4). Thus, the broad prior will contain a fraction of neutron stars rotating impossibly fast for some equations of state. In practice, this issue is not relevant as we find zero posterior support for such rotation rates.
Finally, we also present results where we compare the hypothesis that both bodies are lowmass black holes against the hypothesis that one, or both, bodies are neutron stars, but with a non-specific equation of state model. In this case we are sampling over the tidal deformability of the neutron stars (and deriving all other equation of state specific quantities from this value as explained below). For this analysis we use the same prior choices as above, specifically considering the broad mass prior, except for the following choices on spins and tidal deformability distributions. Here, a black hole is defined as having a spin magnitude uniformly distributed in [0, 0.89] and zero tidal deformability, while a neutron star has a spin magnitude uniformly distributed in [0, 0.05] and a tidal deformability uniformly distributed in [0, 5000]. The upper spin limit for black holes differs from the broad prior considered above. The specific models we study in this case are: where here χ i and Λ i denote the dimensionless spin and tidal deformability of the body i respectively. The analysis of model (iv) allows the two component tidal deformabilities to vary independent of each other. In contrast, the analysis of model (v) uses an equation of state-independent relation between the two tidal deformabilities given the mass ratio of the system [48,49].

Waveform models for the gravitational-wave signal
Gravitational-wave signals from binary neutron stars are distinct from signals from binary black holes. The deformation of the neutron star due to tidal effects and a spin-induced quadrupole, as well as the post-merger behavior of the remnant, will all leave an imprint in the emitted gravitational waves. Unfortunately, no waveform model yet exists that can consistently model all three of these effects, as modeling the post-merger behavior requires expensive numerical simulations. Such simulations have been performed by a number of groups in recent years, enabling an exploration of the expected features of post-merger emission [70][71][72]. Though this work has not yet reached the stage where one can predict waveforms with arbitrary masses, spins and equation of state, unmodeled analyses can target the expected post-merger signal [73][74][75]. Since the post-merger component of the signal occurs at frequencies ∼2000 Hz, where Advanced LIGO and Advanced Virgo sensitivity is decreasing, observing this regime was not possible for GW170817 [3,6,22,76]. For the present analysis, it is therefore safe to ignore the post-merger signal.
We use waveform models that predict gravitational-wave signals for a large range of parameters and include both the effects of the tidal deformation and of the spin-induced quadrupole. In contrast to previous analyses [3,6], where the individual tidal deformabilities Λ of the neutron stars were treated as independent intrinsic parameters, we assume a specific equation of state model and the tidal deformability is computed as a function of the masses. To construct the functions Λ(m), the Tolman-Oppenheimer-Volkoff differential equations [77,78] for the stellar structure are solved numerically together with those for the tidal deformability [12,79].
The quadrupole-monopole effects are caused by the quadrupolar deformations induced by the spin [11], and are quadratic in the spins. These effects are parametrized by a quadrupolemonopole parameter κ. Similarly to the tidal deformability, this parameter is treated as a pure function of the masses κ(m), determined in practice from the value of Λ via fits of approximate universal relations for neutron stars [48].
The waveform models used in this study are essentially the same as in the previous analyses from the LIGO-Virgo collaboration [3,6,39], to which we refer the reader for more details. The TaylorF2 model uses purely post-Newtonian information to generate closed-form waveforms in the Fourier domain, only allowing for spins aligned with the orbital angular momentum. It includes point-particle and spin contributions (see, e.g. [80]) and tidal contributions. In contrast to [6], we consider several variants of TaylorF2, incorporating tidal terms up to the 6, 7 or 7.5 post-Newtonian order [13, 14, 81] 196 .
A closed-form tidal approximant, denoted NRTidal, calibrated to binary neutron star numerical relativity simulations was presented in [82,83]. In this paper, EOB+NRT and Phenom+NRT will denote the waveform models constructed from adding this numer ical relativity tuned approximant for tidal effects onto the underlying point-particle waveform models, SEOBNRv4_ROM [84] and IMRPhenomPv2 [85][86][87][88] respectively. While the EOB+NRT model is restricted to aligned spins, Phenom+NRT incorporates an effective representation of precession effects [85,86]. All three models TaylorF2, Phenom+NRT and EOB+NRT include post-Newtonian quadrupole-monopole terms, up to the third PN order in the phasing [11,89].
In addition, we use two time-domain effective one body models that incorporate tidal effects, SEOBNRv4T [90,91] and TEOBResumS [92]. SEOBNRv4T includes dynamical tides and the effects of the spin-induced quadrupole moment. TEOBResumS incorporates a gravitational-self-force re-summed tidal potential and the spin-induced quadrupole moment. Both models have been shown to be compatible with state-of-the-art numerical simulations of binary neutron stars up to merger [92,93]. Also compare [94,95] regarding recent extensions.
Lacking a prescription for the post-merger signal, all waveforms have an equation of state dependent termination frequency that is close to the merger frequency. For TaylorF2, the termination frequency is determined according to equation (13) of [45] 197 . For EOB+NRT and Phenom+NRT, the termination frequency is computed following equation (11) of [83]. Thus the NRT waveforms with zero tidal deformabilities are not strictly speaking equivalent to the binary black hole waveforms: in the NRT prescription the merger is cut instead of terminating with a ringdown signal as would be the case in the underlying point-particle model. However, this difference only concerns high frequencies and should not affect our evidence computations. For most of the configurations considered in this analysis, a remnant black hole would have a ringdown frequency above ∼ 4 kHz, justifying the above statement.
It is important to note that all our waveform models are imperfect and come with systematic errors. Without a quantitative model to represent this as a modeling uncertainty, we cannot marginalize over the waveform error in our analysis. Comparisons between different waveform models [83,96,97] indicate that TaylorF2 can produce small differences from Phenom+NRT and EOB+NRT in recovering the tidal parameters, while yielding comparable results for the masses and spins. This was indeed observed in the analysis of GW170817 [3,6,38,39], where the TaylorF2 model yielded larger upper limits for the tidal parameters. In order to obtain a rough estimate how differences between waveform models affect the Bayes factor calculations, we compare results obtained with different waveform families.

Equations of state
For this work, we have collected 24 different prescriptions describing neutron star equations of state from the literature, computed by various nuclear physics approximations. The full list together with the relevant citations is given in appendix A.
All equations of state describe zero-temperature matter in β-equilibrium. We require equations of state to respect causality, i.e. the maximum sound speed inside stable neutron stars may not exceed the speed of light at any mass up to the maximum allowed. For compariso n with previous studies [6,39], however, we also include some results for the WFF1 [98] equation of state model, which violates causality before reaching the central density of the maximum-mass neutron star. One equation of state model (HQC18 [9]) incorporates a trans ition between a hadronic phase in the crust and a quark matter phase in the core. The other equation of state models describe purely hadronic matter without significant phase transitions. We deliberately do not use any constraints from previous studies of the detected event [3,6,38,39] to further narrow the equation of state model selection.
The equation of state properties most important for our analysis of the gravitationalwave signal are the maximum gravitational mass of non-rotating neutron stars and their tidal deformability as a function of gravitational mass. In order to derive implications for the merger remnant, we also require the maximum mass of uniformly rotating neutron stars, as well as the maximum rotation rate and moment of inertia as a function of mass. The maximum dimensionless spin for a given mass is also relevant with regard to the spin prior. We compute all those properties for all equations of state used in this work. Key quantities are tabulated in appendix A.
The properties of non-rotating neutron stars which follow from the equations of state used here are shown in figure 1. For comparison, we also show the credible intervals for the component gravitational masses derived in [6] for different priors on the spin. The maximum gravitational mass covers a range 1.92 M (BHF_BBB2) to 2.75M (MS1_PP, MS1B_PP). The discovery of neutron star PSR J0348+0432 [26] with mass 2.01 ± 0.04 M , and the very recent observation of PSR J0740+6620 [99] with mass of 2.14 +0. 10 −0.09 M , provide a lower limit for the maximum mass of neutron stars. However, we also included some models with maximum mass slightly below those constraints in our comparison. Although causality would allow maximum masses even larger than our most massive models (see [100,101]), the chosen range covers most nuclear physics estimates. The smallest tidal deformability within the mass credible interval [6] for the low spin prior is found for the WFF1 equation of state model, with Λ = 61, and the largest one for the MS1_PP equation of state model, with Λ = 3728. In the same mass range, we find that the dimensionless spin at the mass-shedding limit is between 0.66 (KDE0V1) and 0.76 (SKI5).

Gravitational-wave data
For our analysis, we follow [7] and use data from the two LIGO observatories with the full frequency-dependent calibration uncertainties described in [102,103] and following the subtraction of independently measured noise sources as discussed in [104]. As discussed in [3], a non-Gaussian transient occurs in the data recorded by the LIGO Livingston observatory. This non-Gaussian transient is removed from the data as in [3] using the methods discussed in [73]. We have checked that the instrumental transient subtraction does not bias the results of our inference substantially [105]. We also follow [7] in using data from the Virgo observatory using the calibration and noise-subtraction model discussed in [106]. For both the LIGO and Virgo observatories we analyze data between 23 Hz and 2048 Hz.

Evaluating the Bayes factors
In this work we utilize two different strategies to compute the Bayes factors. We use pre-existing tools in the LALInference [42] and RIFT [43] parameter inference codes to compute the evidence for each equation of state model (equation (2)). For LALInference, this involves performing a distinct detailed Bayesian analysis of the data for each proposed model, using nested sampling [107] and thermodynamic integration [108] to compute the evidence [42].
Both methods have been extensively tested [42] and produce consistent results [109].
For RIFT, this calculation is organized as a two-stage process. The first stage estimates a single marginal likelihood P(s|x, Λ 1 , Λ 2 , I) = dθP(s|θ, x, Λ 1 , Λ 2 ; M i , I)P(θ) as a function of the extrinsic parameters θ and the intrinsic binary parameters x = (m 1 , m 2 , χ 1 , χ 2 ) and (Λ 1 , Λ 2 ) without assuming an equation of state. The next stage computes the evidence via where Λ(m|M i ) denotes the deformability as function of mass assuming the equation of state model M i . Using a precomputed marginal likelihood, this second stage can be performed very rapidly for any equation of state, in seconds to minutes depending on the accuracy required. Unless otherwise noted, all RIFT calculations are performed using the marginalized likelihoods P(s|x, Λ 1 , Λ 2 ) employed in [7].

Results
For each equation of state described in section 2.4, the Bayes factors are computed using either LALInference or RIFT, using both the narrow and broad priors discussed in section 2.6. The Bayes factors for all waveform approximants and the statistical uncertainties are given in appendix B. The Bayes factors calculated adopting the narrow prior are shown in  proper circumferential radius. All sequences terminate at the maximum mass model (except WFF1, which is cut where the maximum sound speed reaches the speed of light). The legend entries are sorted by the radius of a 1.36 M neutron star. For comparison, we show the 90% symmetric credible interval given in [6] for the component masses, with dashed vertical lines for the low-spin prior, dotted lines for the high spin prior, and the solid line for the equal-mass limit 1.36 M . The shaded vertical area marks the measured mass for pulsar PSR J0348+0432 [26]. tables B1 and B3 for LALInference and RIFT, respectively. Similarly, for the broad prior, the LALInference results are given in table B2 and the RIFT results in table B4. The Bayes factors for the LALInference runs are quoted with respect to the low-mass binary black hole hypothesis. The RIFT results re-use the analysis presented in [7] and are not normalized with respect to the BBH case due to the way in which the fits are performed. Instead, the RIFT Bayes factors are quoted relative to the arbitrary hypothesis that GW170817 was a binary neutron star with the SLY9 equation of state, which has modest support and shows relatively little dispersion between waveform approximants in the LALInference results.
For the narrow prior, the gravitational-wave data does not favor and nor can it rule out any single equation of state with respect to the low-mass binary black hole hypothesis. This can be seen in figure 2, although the data clearly favors softer equations of state, which predict more compact neutron stars with lower tidal deformabilities. The broad prior case, shown in figure 3, exhibits the same trend. The relation is visualized more directly in figure 4, showing fiducial neutron star properties as a function of the Bayes factor with respect to the KDE0V model, which is the most favored for the narrow prior with TaylorF2 (7.5pN) approximant. The models with largest fiducial radius and tidal deformability are very strongly disfavored (MS1_PP) or at least strongly disfavored (MS1B_PP) compared to the numerous models with small fiducial radius and deformability, such as KDE0V. Figure 4 also shows upper limits for radius and tidal deformability given by parameter estimation studies [6,38,39]. For all our equations of state, and mass ratios in the range q > 0.7, the effective tidal deformability and the tidal deformability at the fiducial mass agree better than 6%. To this accuracy, we can compare the upper limits to the values plotted for the fiducial mass models. We find that some of our models with tidal deformabilities above the limits are not strongly disfavored with respect to the most favored equations of state.
Although we can only compute Bayes factors between the models at hand, the Bayes factors exhibit a trend of decreasing with an increasing fiducial radius and tidal deformability. If this trend holds in general, then for a fiducial mass of 1.36M models with tidal deformability Λ(M = 1.36) 1600 or radius R(M = 1.36) 14.9 km are very strongly disfavored. For the broad prior, the gravitational-wave data alone does not prefer any single equation of state, as was the case for the narrow prior. However, MS1_PP and MS1B_PP are now strongly disfavored against the low-mass binary black hole hypothesis, as shown in figure 3. Normalizing the Bayes factors against KDE0V, we find that MS1_PP is very strongly and MS1B_PP strongly disfavored with respect to the softer KDE0V equation of state, as shown in figure 4.
In contrast to the narrow prior, we find that for some equations of state there is significant posterior support for a mixed binary where the heavier object is a light black hole. The reason for this is that the mass range is less constrained and exceeds the maximum possible mass of non-rotating neutron stars for those equations of state. Models BHF_BBB2, KDE0V1, HQC18, and KDE0V have 13, 7, 5, and 7 percent posterior probability, respectively, for a mixed binary. The remaining models have support below 5% for mixed binary, in particular MS1_PP and MS1B_PP have zero support for mixed binaries. For the broad prior and the WFF1 equation of state, there is a ≈14% posterior probability for a central neutron star density in an unphysical (causality violating) region of the parameter space.
For both sets of priors, we find that the bulk of equations of state considered in this paper produce evidences that are consistent to within an order of magnitude. This is in line with the expectations for second-generation gravitational-wave detectors, where ∼O(20) detections will allow hypothesis ranking to distinguish between hard, moderate and soft equations of state [44]. Even if the true equation of state were to be one of those considered in this paper, we should not necessarily expect this to be the most preferred equation of state in the hypothesis ranking [45]. For example, this could arise due to the effects of noise or other systematics. More realistically, it is also possible that the true equation of state is not contained in the finite set of models considered here, though we expect the highest ranked hypotheses to share similar macroscopic features to the true neutron star equation of state.
It is important to assess the systematic bias inherent in our waveform models. Such systematic biases could arise from a different treatment of the underlying point particle limit, which can start to become significant at higher frequencies where tidal effects also become more prominent. Additionally, the different waveform models employ different prescriptions for the tidal contribution to the phasing, which can lead to discrepancies in the recovery of tidal parameters. As demonstrated in previous studies [39,83,96,97,110], the tidal effect calculated in the TaylorF2 waveform models is typically smaller than that predicted using the  NRTidal prescription. As such, the tidal parameters inferred from the data using the TaylorF2 waveform model will typically be larger in order to compensate for the reduced tidal effect. Consequentially, we find that the Bayes factors calculated using the TaylorF2 model tend to show greater support for stiffer equations of state, which yield larger tidal deformabilities, in comparison to the effective-one-body or phenomenological waveform models. However, as can be seen in figures 2-4, the Bayes factors between the different waveform approximants still show reasonably good agreement. Finally, we also ask if can we rule out the binary black hole hypothesis by computing Bayes factors for the binary neutron star hypothesis, or the hypothesis that there was only one neutron star in the binary, without assuming a specific equation of state. We consider both for both priors and all waveform approximants. For comparison, we plot limits from the literature. The solid green horizontal lines mark the upper limit for the coalescing neutron star radii from [39]. Red dashed lines mark the upper limit given for the radius given in [38] (assuming a common radius for the two stars) and the one-sided upper limit for the effective tidal deformability from [38]. The blue dash-dotted lines mark the upper limit (highest posterior density interval) for the effective tidal deformability from [6]. the hypotheses that the single neutron star is heavier and lighter than its potential black hole companion. The specific priors that we use in these cases were discussed earlier in section 2.2. Table 1 presents the Bayes factors compared to the 'binary neutron star with independent equation of state' hypothesis. We present results with different waveform approximants and methodologies to compute the evidence of each model.
From the results in table 1 we can draw a number of conclusions. First, we find that the computed Bayes factors agree closely between different waveform approximants, indicating that uncertainties in the waveform modeling are not compromising our results. Second, we find that both nested sampling and thermodynamic integration return consistent results for the Bayes factors, giving us confidence on the validity of the calculation. Third, and most interesting from a physical perspective, we find that the GW170817 data alone does not strongly favor one hypothesis over the other. All of the Bayes factors we compute are in the range of 1 to 5, suggesting that the gravitational-wave data alone do not show a strong preference for any of the hypotheses considered here.
It may seem strange that the binary neutron star model, assuming no common equation of state, is the least favored hypothesis. Bayes factors depend on the Occam penalty, i.e. how well a model fits the given data (goodness-of-fit) and how many constrained parameters a model includes (dimensionality). We expect all models to be preferred over the binary neutron star model, assuming a common equation of state, if they fit the data equally well. The binary black hole model has two fewer parameters and the neutron star-black hole hypotheses have one fewer parameter. A brief back-of-the-envelope calculation can be used to support the Bayes factors seen in table 1. The Occam penalty is given by δθ/∆θ, the posterior width over the prior width, and penalizes a model for each additional constrained parameter it employs, since δθ < ∆θ. From section 2.2, we see that the width of the prior on the tidal parameters is ∆Λ ∼ 5000 whereas the width on the posterior is on the order of δΛ ∼ 1000 [6]. From this, we expect the binary neutron star hypothesis with a common equation of state to be preferred over an independent equation of state hypothesis by a factor ∼5, as it uses one fewer parameter. This is in agreement with the results found in table 1 and as demonstrated in [49] for simulated binary neutron star signals.

Implications
After comparing the relative probabilities of the different equation of state models, we now discuss the implications that follow from assuming any of the equations of state are correct. We also discuss constraints for the maximum mass of neutron stars that follow from assumptions on the fate of the remnant, but without assuming any equation of state.

Baryonic mass and central density
The total baryonic mass (defined here as baryon number times a fiducial mass constant 1.66 × 10 −27 kg) is the most important quantity for the evolution during the post-merger phase, including the lifetime of the remnant and the mass of the accretion disk, which are relevant for the short gamma-ray burst and the kilonova. The maximum density inside the initial neutron stars is important to interpret our results properly. The relative probabilities we compute for different equation of state models only signify that in the density range occurring in the coalescing neutron stars, matter is more likely described by one equation of state model than another.
To compute the baryonic mass, we start with the posterior samples obtained assuming a given equation of state model. From the gravitational masses and spins in a given sample, we obtain the baryonic masses of the neutron stars. Those are computed using interpolation tables we build for sequences of non-rotating and maximally rotating neutron stars for each equation of state model. To account for the spin, we expand the gravitational mass at given baryonic mass to second order in the squared dimensionless spin parameter. The zerothand first-order coefficients are determined from the binding energy and moment of inertia of the non-rotating model. The second-order coefficient is then determined from the angular momentum and gravitational mass of the maximally rotating model. The spin corrections are relatively small compared to the statistical error both for narrow and broad priors. To obtain limits for the energy densities that can be reached before merger, we use the central densities of non-rotating neutron stars with baryonic masses computed as described above.
The results for the narrow prior are given in table 2. We find that the pre-merger gravitational-wave signal does not carry information about matter at energy densities above 0.5-1.3 × 10 15 g cm −3 , assuming the true equation of state is similar to one of the models considered here. This agrees well with the results by [39] (figure 2 therein) obtained using parametrized equations of state, and also with results in [40] based on Gaussian process models for the equation of state.

Fate of merger remnant
Next, we assess the fate of the merger remnant, which could either form a black hole directly at merger (prompt collapse), a short-lived hypermassive neutron star, a long-lived supramassive neutron star, or even an unconditionally stable neutron star. This is particularly relevant with regard to the observed short gamma-ray burst and kilonova counterparts. One possible engine for short gamma-ray bursts is given by a black hole submerged in a massive disk [52]. Numerical simulations suggest that prompt black hole formation at merger might be in tension with the production of a short gamma-ray bursts [29,112], favoring delayed collapse. The other model for the engine consists of a rapidly rotating and strongly magnetized neutron star (magnetar scenario [53,113,114]) embedded in a disk. In this model, the remnant needs to survive at least long enough to initiate a jet (however, there can be a delay until the jet produces the observed short gamma-ray burst). The fate of the remnant is also tightly related to the ejection of matter and the spectral evolution of the resulting kilonova [5,[115][116][117][118]. For example, the amount, composition and velocity of ejected matter required for the observed kilonova might be in tension with expectations for the prompt collapse case (for a recent review on multimessenger constraints, see [119]). However, the theoretical modeling of the kilonova is a very complex and rapidly evolving field of research (e.g. [120][121][122][123] and the references therein). Table 2. Implications derived from the equation of state model specific parameter estimation runs with narrow prior and TaylorF2 (7.5 PN) waveform model. M 0 B is the total initial baryonic mass. We provide the median value and the 90% two-sided credible bounds. E c is the 90% one-sided upper limit for the central energy density of the heavier neutron star. The normalized mass m 0 is simply a linear transform of M 0 B defined by equation (4), such that the supramassive mass range for the given equation of state model corresponds to m ∈ (0, 1). The potential impact of mass ejection is given in terms of ∆ 0.1 =m 0 −m R , computed from equation (5) for a fiducial ejecta mass 0.1 M (the given value is the 90% one-sided upper limit). The 'prompt collapse' column denotes whether a black hole can be formed directly during the merger, based on available estimates of the corresponding mass threshold, and assuming prompt collapse only occurs for hypermassive systems. Hypermassive cases lacking data on the prompt collapse threshold are denoted with a dash, while 'possibly' refers to the case where the available bounds on the threshold are compatible with both outcomes. I max and F max are upper limits (see main text) for moment of inertia and rotation rate for a uniformly rotating remnant. I max is given in geometric units of 4.335 871 × 10 43 g cm 2 . The most relevant quantity regarding the remnant type is the baryonic mass M R B of the system, in comparison to the maximum baryonic mass of non-rotating neutron stars, M TOV B , and the maximum baryonic mass of neutron stars uniformly rotating at the mass-shedding limit, M ROT B . Since the baryon number is conserved, the remnant mass is given by The quantity m R serves as a rough indication for the remnant lifetime τ R . A stable remnant will never form a black hole, a hypermassive remnant either forms a black hole directly at merger or collapses within a few hundred ms. Within the narrow supramassive mass range, the lifetime decreases sharply. Computing τ R as a function of m R , however, is a difficult problem beyond the scope of this paper. Although we expect m R to be the most important parameter, note that τ R might also depend on other factors such as mass ratio, spins, and equation of state.

Equation of state
In table 2, we provide the value m 0 corresponding to the median value of the total initial baryonic mass inferred for the narrow prior. Applying equation (4) In [116], the ejecta mass was estimated from the observed kilonova to around 0.078 M . If we assume more conservatively that M E B < 0.1 M , we find that 0 <m 0 −m R < 0.3 for the equations of state considered here (see table 2). The above ejecta mass limit should be regarded as a fiducial value since modeling the ejecta mass, either by numerical simulations or inferred from the kilonova, is an ongoing field of research.
For the narrow prior, we find that the MS1_PP, MS1B_PP, and MPA1 equation of state models lead to a stable remnant and are therefore incompatible with the black hole with disk short gamma-ray burst scenario (technically, the MPA1 also allows marginally supramassive remnants, which can for all practical purposes be considered stable). The H4, SLY, KDE0V, KDE0V1, SKOP, SLY2, BHF_BBB2, and HQC18 equation of state models result in a shortlived hypermassive remnant. For APR4_EPP, RS, SLY9, SLY230A, SKI2, SKI4, SK255, and SKMP, the remnant is in the upper supramassive or lower hypermassive mass range, and we cannot predict if the lifetime is longer or shorter than the observed short gamma-ray burst delay. For SKI3, SKI5, SKI6, and SK272 the remnant mass is in the supramassive range with m R < 0.9, and likely long-lived.
We now turn to discuss whether the merger can lead to prompt black hole formation for a given equation of state model, based on gravitational-wave data alone. Prompt collapse occurs above some critical mass, which depends on the equation of state. Known examples of prompt collapse in numerical binary neutron star simulations generally involve hypermassive systems. If a parameter estimation run for a given equation of state model yields posterior support only in the supramassive and stable mass range, we assume that no prompt collapse can occur. For hypermassive systems, we use existing thresholds [18] from numerical simulations, given in terms of total gravitational mass of the binary. We consider a systematic error of 0.1 M due to the granularity of the simulated systems. The available thresholds are for equal-mass, nonspinning systems. To account for a possible impact of mass ratio and spin, we further increase the error by 0.05 M . This is a conservative estimate assuming that the main influence is by the change in dynamically ejected mass, which has been computed for many models in [124]. Finally, we express the thresholds for equal mass systems in terms of total baryonic mass and use it as an approximate threshold on M tot B for the generic case. For the narrow prior, we can rule out prompt collapse for the APR4_EPP, H4, MPA1, MS1_PP, MS1B_PP, SKI3, SKI4, SKI5, SKI6, and SKI272 equation of state models. The SLY equation of state model is compatible both with prompt and delayed collapse (due to the systematic error of the threshold mass from numerical relativity). For the remaining equations of state, no data on prompt collapse threshold is available.
With improved theoretical modeling, the wealth of electromagnetic observations may yield reliable constraints on the lifetime of the remnant. If prompt collapse could be ruled out, the values of M 0 B obtained from gravitational-wave data could be taken as a lower limit for the prompt collapse threshold. Considering a potential impact of the mass ratio on the actual threshold, one has to take into account that the limits on M 0 B are the result of a marginalization over mass ratio.

Constraining the maximum neutron star mass
Within the black hole with disk scenario of short gamma-ray bursts, the remnant cannot be stable. This was already exploited in [4] to derive limits on the maximum mass of non-rotating neutron stars. The delay of 1.7 s [4] between merger and short gamma-ray burst would also constrain the lifetime of the remnant. The relatively short lifetime indicates that the remnant either cannot be supported by uniform rotation alone (hypermassive remnant), or that it is at least close to the critical mass (how close exactly is an important open question).
In the following, we estimate the maximum mass of non-rotating neutron stars under the assumption that a black hole was produced in the binary neutron star merger before the observed short gamma-ray burst. For this, we first rewrite equation (4) aŝ From table 2, we find that M 0 B varies only slightly with equation of state. Therefore, we take the bound M 0 B < 3.05 M for the equations of state considered here as a generic bound. In addition, we assume k m < 7, as fulfilled by all equations of state considered here. The remaining open question is below which value m R the remnant lifetime is above the short gamma-ray burst delay. We already know that the lifetime is infinite unless m R > 0. Future studies might provide tighter bounds, which will tighten the upper bound obtained for M TOV B . For example, if we demand m R > 0.4, then M TOV B < 2.9 M . Another potential method to limit the maximum mass of non-rotating neutron stars is to limit the angular momentum that needs to be removed before collapse can occur, and which becomes zero for m R = 1. In [27], it was argued that the remnant mass must be close to the hypermassive range since the large spindown energy that would otherwise have been deposited in the environment would be incompatible with the observed counterparts (however, there are observational gaps [5], e.g. in x-rays, which could lead to energy output being missed). This assumption corresponds to m R 1. From equation (6), we would obtain a corre sponding limit M TOV B < 2.67 M (for the narrow prior). If we further assume that the maximum baryonic mass is 15%-23% larger than the maximum gravitational mass, as fulfilled for all equation of state models considered here, we obtain a limit M TOV G 2.32 M . Based on additional assumptions, excluding prompt collapse and using an approximate expression for the prompt collapse threshold ignoring systematic errors, [27] provides a limit on the maximum neutron star mass which is 6% tighter. Reliable prompt collapse thresholds for a representative selection of equation of state models might hence be useful in this respect.
In [29], a limit M TOV G 2.28 ± 0.23 M was provided, again based on the assumption of having a hypermassive remnant. The differences to our analysis are mainly the use of gravitational masses instead of conserved baryonic masses, and a more constraining assumption regarding the range of ratios between maximum gravitational masses for uniformly rotating and non-rotating neutron stars. Such effects are included in the systematic error given in [29] for their maximum mass limit, which is therefore compatible with our result.
A similar calculation can be found in [54], yielding a limit M TOV G 2.15-2.25 M . It is based again on the detected limits for the total initial gravitational mass, and further assumes a total energy loss of (0.15 ± 0.03) M , as well as a fixed difference 0.4 M between maximum gravitational masses of uniformly rotating and non-rotating neutron stars. These assumptions might explain the tighter limit compared to our result (for example, we find M ROT G −M TOV G = 0.34 M for the BHF_BBB2 equation of state). A recent analysis [125], which is not relying on assumptions about m, indicates that the remnant at the time of collapse might not rotate rapidly, and provides a limit M TOV G 2.3 M . By using the detected gravitational-wave signal alone, without assumptions on the remnant fate, [40] inferred a limit M TOV G 2.46, compatible with our results. This is based on modeling the equation of state using Gaussian processes, trained on a set of seven equations of state and respecting causality by construction. The Gaussian process model is further informed by the gravitational-wave data and then used to extrapolate the equation of state from densities relevant for the inspiral signal up to those required to determine the maximum mass.

Impact of prior and waveform model
In order to assess the waveform systematics, we repeat the computation of the remnant properties for a representative equation of state model using all waveform approximants. The results are shown in table 3. The differences for M B and m 0 are small compared to the statistical uncertainties, and negligible with regard to our conclusions on remnant classification, prompt collapse, and the maximum neutron star mass. This is to be expected, since those results depend directly only on the initial neutron star masses, but not on the tidal deformabilities, which are most sensitive to the waveform model. Table 4 shows the results obtained using the broad prior restricted to the binary neutron star case. The restriction is necessary as the heavier object exceeds the maximum mass allowed for a neutron star for parts of the broad prior. Comparing to the narrow prior in table 2, we find that the total initial baryonic mass distributions are broadened towards higher masses. This is a consequence of the wider posterior distribution for the mass ratio; for fixed chirp mass, the total baryonic mass for unequal mass binaries is higher than for equal masses.
The broader distribution leads to differences regarding the fate of the remnant. While we ruled out prompt collapse for the H4 and APR4_EPP equations of state with the narrow prior, the respective total mass posteriors for the broad prior also have support in the mass range where prompt collapse is allowed according to the thresholds from numerical relativity. This result might still change with more accurate numerical-relativity results for the threshold, given its large systematic error.
The KDE0V, SKI3, SKI4, SKI5, SKI6, and SK272 equations of state, which lead to supramassive remnants assuming the narrow prior, might also lead to hypermassive remants when assuming the broad prior. The MS1_PP, MS1B_PP equations of state still lead to stable remnants, while the MPA1 equation of state model can also result in a long-lived supramassive remnant for the broad prior. The bound on the maximum Tolman-Oppenheimer-Volkoff mass when assuming the black hole with disk short gamma-ray burst model becomes less strict as well. Using only m R 0, we obtain M TOV B < 3.4M . Assuming m R 1, we would obtain M TOV B < 2.98M . Another consequence of the larger mass ratios is that the mass posterior of the heavier neutron star extends to larger values. Therefore, higher central energy densities are reached. Depending on the equation of state, the inspiral gravitational-wave signal might be sensitive to the equation of state at energy densities up to 2 × 10 15 g cm −3 (compare table 4).

Remnant rotation rate
In case the remnant is long-lived, the maximum rotation rate and moment of inertia are of interest for searches of long-lived gravitational-wave signals from the remnant (see [20,22]), and might also be relevant for modeling the magnetar short gamma-ray burst scenario. In the following, we estimate upper limits for rotation rate and moment of inertia.
The total angular momentum at merger is not known, although numerical simulations indicate that it typically exceeds the maximum amount that could be realized in a uniformly rotating system of the same mass (see [21]). However, the disk can account for a significant fraction of the total angular momentum. Most of the disk angular momentum is transported outwards during accretion (together with a fraction of the mass). The angular momentum of the final remnant could depend on the details of the accretion process and might be less than the critical one (the opposite assumption is made in [21]), and it can decrease further on longer timescales, e.g. by magnetic spindown. The rotation rate of models near the maximum mass can increase with decreasing angular momentum, because the compactness also increases. The maximum rotation rate can in principle exceed the rotation rate at the mass-shedding limit. Here, we use the rotation rate at the mass-shedding limit as an estimate for the maximum rotation rate of the remnant.
From the mass posteriors for each equation of state model, we compute the posteriors for the rotation rate and moment of inertia of uniformly rotating neutron stars with baryonic mass equal to the total baryonic mass and maximum angular momentum. Tables 2 and 4 show the numerical results for the narrow and broad priors. For both, we find a maximum rotation rate below 1.9 kHz for all equations of state. We also provide values for remnants in the hypermassive mass range. For those cases, we limit the baryonic mass above to the maximum allowed for uniform rotation. The numbers we obtain can be regarded as estimates for the hypothetical case that the excess mass was somehow lost during merger, or remains in a disk.

Conclusions
In this work, we employ Bayesian model selection as a tool to discriminate between different theoretical models for the neutron star equation of state using only the gravitational-wave data. This allows us to make quantitative statements about the relative likelihood of a given equation of state model. This work complements the analysis of tidal effects of GW170817 in [6,38] and the analysis of the neutron star radii and equation of state in [39]. As in previous studies [6,38], we assume a priori that the source of GW170817 resides in NGC 4993 and adopt fixed sky location consistent with electromagnetic observations. The distance to GW170817 is determined from the gravitational-wave data alone, as in [6]. We do not allow tidal parameters to vary independently and instead enforce that each neutron star obeys a given equation of state model. The Bayesian evidence is computed using two independent methods that show good agreement and demonstrate the robustness of our results. Under minimal assumptions on the astrophysical priors for the component masses and spins, we find that the models predicting the largest fiducial radii and tidal deformabilities are very strongly disfavored (MS1_PP) or at least strongly disfavored (MS1B_PP) with respect to the equations of state that predict a small fiducial radius and tidal deformability, such as KDE0V. In addition, we find that for both priors considered in this work the gravitational-wave data does not favor and nor can it rule out any particular equation of state against the low-mass binary black hole hypothesis.
In order to gauge waveform systematics, we compare the Bayes factors using different waveform approximants. The results show good agreement between the different waveform models, providing confidence that systematic uncertainties are small compared to statistical uncertainties.
In a broader context, we also address the question as to whether we can exclude the possibility of GW170817 being a binary black hole by computing Bayes factors against the binary neutron star and neutron star-black hole hypotheses. From the gravitational-wave data alone, no single hypothesis is strongly preferred and we are unable to exclude the binary black hole hypothesis. We find that the Bayes factors computed using different waveform approximants are in good agreement, supporting the notion that waveform systematics are small compared to statistical errors.
For some equation of state models, we can predict the outcome of the merger; for two models the remnant has to be a stable neutron star, for some it will be a long-lived supramassive neutron star, and for some a short-lived hypermassive neutron star or black hole. We can further rule out prompt collapse for a number of models. Considering the Bayes factor for each model, we find that only indefinitely stable remnants are disfavored. Based on gravitationalwave data alone, all scenarios from very long-lived remnants to prompt black hole formation at merger remain possible. Further constraints on the outcome will likely have to rely on the electromagnetic counterparts and modeling of the kilonova and the short gamma-ray burst. We find that making the assumption of a short-lived remnant-as required by some short gammaray burst models-would allow us to put further constraints on the equation of state model. In particular, we improve the upper limit on the maximum mass of non-rotating neutron stars, compared to an earlier study [4] (which is based on an analysis that did not assume that both neutron stars obey the same equation of state). We also find that three equation of state models are ruled out if a black hole was formed before the short gamma-ray burst, including the two disfavored by gravitational-wave data alone.
Research, for the construction and operation of the Virgo detector and the creation and support of the EGO consortium. The authors also gratefully acknowledge research support from these agencies as well as by the Council of Scientific

Appendix A. Equations of state
In the following, we provide details for the equations of state employed in this work, all of which represent estimates from nuclear physics. They employ different estimates for the forces between constituents, for example from nucleon-nucleon scattering experiments. Those forces are often parametrized (e.g. Skyrme forces) before use in equation of state computations. Other important observable parameters to calibrate equation of state models are the nuclear symmetry energy and incompressibility. Further, the many-particle quantum field theory problem is approximately solved using different schemes, such as liquid-drop model, relativistic/nonrelativistic mean field (MF) theory, Hartree-Fock (HF) method, Brueckner-Bethe-Goldstone (BG) methods, or Wigner-Seitz (WS) cells. Some methods apply the simplification of a single microstate, while others consider nuclear statistic equilibrium. One important aspect is which types of constituents are considered, for example if hyperons are included or not. Most equation of state models considered describe hadronic matter, except the HQC18 equation of state model, which assumes a hadronic crust around a quark matter core [9].
For details on the various models, we refer to the references given in table A1, which list all our equation of state models, the respective sources for the numeric tables, and the original publications. Further, we provide the maximum gravitational and baryonic masses of non-rotating neutron stars, and the maximum baryonic mass of uniformly rotating neutron stars. The latter were obtained using the Rapidly Rotating Neutron Star (RNS) code described in [126].
Most of our equations of state were obtained from three catalogs. The CompOSE project [127] provides nuclear physics equation of state data and a supporting software framework.
We used this public code (selecting the first-order interpolation option and a resolution of 500 points) to create tables. Another collection of tables is provided by [128]. To obtain the final, continuous equations of state, we employ two different interpolation schemes. The results in section 4 and table A1 are computed using monotonic spline interpolation (PCHIP), while the equation of state related code used during parameter estimation employs standard cubic spline interpolation instead. For all models, we compared sequences of non-rotating neutron stars between the two. Since the low resolution and gaps present in the MS1 and MS1B tables cause significant ambiguities in the tidal deformability and radius, we use well defined piecewise polytropic approximations instead, taken from a third catalog [47] providing analytic fits, and denoted here by MS1_PP and MS1B_PP, respectively. When comparing with [6], note that their WFF1, H4, MPA1, SLY and APR4 models also refer to the piecewise polytropic variants from [47], in contrast to this work.
We applied the following adjustments: in general, we removed isolated data points from the original tables that violate thermodynamic constraints for zero temperature equations of state. For the H4 equation of state model from [128], the low-density part (<10 13 g cm −3 ) was  missing and filled in using the corresponding values of the H3 equation of state model from the same source. For the equations of state from [129], the sound speed becomes superluminal shortly before the maximum mass model is reached, both for the tabulated versions [128] (AP4) or [127] (APR), and the piecewise polytropic fit from [47] (APR4). We therefore use the modified APR4 (denoted APR4_EPP here) from [130], which adds additional segments at high density to ensure causality.

Appendix B. Bayes factors
Here we present the explicit values for the Bayes factors computed using the nested sampling algorithm implemented in LALInference and RIFT. We show results for all of the equations of state discussed in section 2.4 for both the wide and narrow priors outlined in section 2.2.