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

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–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

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle O_{12} = \frac{P(M_1 | s , I)}{P(M_2 | s , I)} = \frac{P(s | M_1, I)}{P(s | M_2 , I)} \frac{P(M_1 | I)} {P(M_2 | I)} = B_{12} \frac{P(M_1 | I)} {P(M_2 | I)} , \nonumber \end{align} \tag{ 1 }$$

where $\mathcal{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₁₂ reduces to the Bayes factor B₁₂. The evidence is calculated by marginalizing over all parameters weighted by their prior probabilities $P(\xi | M_i , I)$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:Z} \mathcal{Z} = P(s | M_i , I) \propto \int P(s | \xi , M_i, I) \, P(\xi | M_i , I) \; {\rm d} \xi , \nonumber \end{align} \tag{ 2 }$$

where the likelihood of obtaining the data realization s given a gravitational-wave signal with physical parameters $\xi$ is

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle P(s | \xi , M_i , I) \propto {\rm e}^ {- \left\langle s - h(t,\xi) \, | \, s - h(t,\xi) \right\rangle / 2}. \nonumber \end{align} \tag{ 3 }$$

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 $\xi$.

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(\xi | 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 $\phi_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–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.33 M_{\odot}$ and a standard deviation of $0.09 M_{\odot}$ [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.7 M_{\odot}$ and $3.0 M_{\odot}$ with the constraint that $m_2 \leqslant m_1$. Priors on the chirp mass $\mathcal{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 fastest-spinning known neutron star, which has a dimensionless spin $\leqslant 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 low-mass 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:

• (i)  
  Binary Black Hole: $\chi_1 \in [0,0.89]$, $\chi_2 \in [0,0.89]$, $\Lambda_1=\Lambda_2=0$.
• (ii)  
  Black Hole—Neutron Star: $\chi_1 \in [0,0.89]$, $\chi_2 \in [0,0.05]$, $\Lambda_1=0$, $\Lambda_2 \in [0,5000]$, $m_{\rm BH} > m_{\rm NS}$.
• (iii)  
  Neutron Star—Black Hole: $\chi_1 \in [0,0.05]$, $\chi_2 \in [0,0.89]$, $\Lambda_1 \in [0,5000]$, $\Lambda_2 =0$, $m_{\rm NS} > m_{\rm BH}$.
• (iv)  
  Binary Neutron Star, no assumptions on equation of state: $\chi_1 \in [0,0.05]$, $\chi_2 \in [0,0.05]$, $\Lambda_1 \in [0,5000]$, $\Lambda_2 \in [0,5000]$.
• (v)  
  Binary Neutron Star, neutron stars obey same equation of state: $\chi_1 \in [0,0.05]$, $\chi_2 \in [0,0.05]$, $\newcommand{\e}{{\rm e}} \Lambda_s \equiv (\Lambda_2+\Lambda_1)/2 \in [0,5000]$, $\newcommand{\e}{{\rm e}} \Lambda_a \equiv (\Lambda_2-\Lambda_1)/2 = \Lambda_a(\Lambda_s, m_2/m_1)$.

where here $\chi_{i}$ and $\Lambda_{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].

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–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–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 $\Lambda$ 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 $\Lambda(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 quadrupole-monopole parameter $\kappa$. Similarly to the tidal deformability, this parameter is treated as a pure function of the masses $\kappa(m)$, determined in practice from the value of $\Lambda$ 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]¹⁹⁶.

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 numerical relativity tuned approximant for tidal effects onto the underlying point-particle waveform models, SEOBNRv4_ROM [84] and IMRPhenomPv2 [85–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]¹⁹⁷. 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 ${\sim}$$4\, {\rm 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.

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 $\beta$-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 comparison 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 transition 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 gravitational-wave 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_\odot$ (BHF_BBB2) to $2.75 M_\odot$ (MS1_PP, MS1B_PP). The discovery of neutron star PSR J0348+0432 [26] with mass $2.01 \pm 0.04 \, M_\odot$, and the very recent observation of PSR J0740+6620 [99] with mass of $2.14^{+0.10}_{-0.09} \, M_\odot$, 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 $\Lambda = 61$, and the largest one for the MS1_PP equation of state model, with $\Lambda= 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).

Zoom In Zoom Out Reset image size

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.

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,\Lambda_1,\Lambda_2, I) = \int {\rm d}\theta P(s|\theta,x,\Lambda_1,\Lambda_2;M_i,I)P(\theta)$ as a function of the extrinsic parameters $\theta$ and the intrinsic binary parameters $x=(m_1,m_2,\chi_1,\chi_2)$ and $(\Lambda_1,\Lambda_2)$ without assuming an equation of state. The next stage computes the evidence via $P(s | M_i , I) = \int {\rm d}x P(x) P(s|x, \Lambda(m_1|M_i),\Lambda(m_2|M_i);M_i,I)$, where $\Lambda(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,\Lambda_1,\Lambda_2)$ employed in [7].

The total baryonic mass (defined here as baryon number times a fiducial mass constant $1.66\times 10^{-27}\, \mathrm{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 zeroth- and 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 \times 10^{15}~\mathrm{g}~\mathrm{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.

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_\mathrm{B}^0$ 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 $\tilde{m}_0$ is simply a linear transform of $M_\mathrm{B}^0$ defined by equation (4), such that the supramassive mass range for the given equation of state model corresponds to $\tilde{m} \in (0,1)$. The potential impact of mass ejection is given in terms of $\Delta_{0.1} = \tilde{m}_0 - \tilde{m}_\mathrm{R}$, computed from equation (5) for a fiducial ejecta mass $0.1\,M_\odot$ (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_\mathrm{max}$ and $F_\mathrm{max}$ are upper limits (see main text) for moment of inertia and rotation rate for a uniformly rotating remnant. $I_\mathrm{max}$ is given in geometric units of $4.335\,871 \times 10^{43} ~\mathrm{g}\, ~\mathrm{cm}^2$.

Equation of state	$M_\mathrm{B}^0$ $(M_\odot)$	$E_c / c^2$ $(10^{15}~ \mathrm{g}~\mathrm{cm}^{-3})$	$\tilde{m}_0$	$\Delta_{0.1}$	Prompt collapse	$I_\mathrm{max}$ $(M_\odot^3)$	$F_\mathrm{max}$ $(\mathrm{kHz})$
APR4_EPP	$3.021^{{+}0.019}_{{-}0.009}$	1.05	$0.883^{{+}0.041}_{{-}0.019}$	0.22	No	84	1.79
MS1_PP	$2.940^{{+}0.014}_{{-}0.008}$	0.47	$-0.608^{{+}0.023}_{{-}0.013}$	0.17	No	153	1.00
MS1B_PP	$2.949^{{+}0.014}_{{-}0.008}$	0.49	$-0.581^{{+}0.023}_{{-}0.013}$	0.16	No	152	1.03
KDE0V	$3.013^{{+}0.019}_{{-}0.009}$	1.17	$1.859^{{+}0.050}_{{-}0.025}$	0.27	—	66	1.86
SKI6	$2.994^{{+}0.016}_{{-}0.009}$	0.82	$0.848^{{+}0.036}_{{-}0.018}$	0.22	No	97	1.53
H4	$2.978^{{+}0.015}_{{-}0.009}$	0.65	$1.385^{{+}0.035}_{{-}0.017}$	0.22	No	104	1.36
MPA1	$3.004^{{+}0.016}_{{-}0.009}$	0.74	$-0.016^{{+}0.029}_{{-}0.016}$	0.18	No	114	1.34
SK255	$2.968^{{+}0.015}_{{-}0.008}$	0.80	$1.120^{{+}0.038}_{{-}0.020}$	0.24	—	91	1.61
SKI2	$2.968^{{+}0.016}_{{-}0.008}$	0.72	$1.025^{{+}0.036}_{{-}0.019}$	0.23	—	100	1.53
SKI4	$2.997^{{+}0.017}_{{-}0.008}$	0.83	$0.911^{{+}0.035}_{{-}0.018}$	0.22	No	94	1.58
SKOP	$2.994^{{+}0.016}_{{-}0.009}$	1.02	$1.808^{{+}0.043}_{{-}0.022}$	0.26	—	73	1.70
SLY2	$3.005^{{+}0.018}_{{-}0.009}$	1.02	$1.401^{{+}0.043}_{{-}0.022}$	0.25	—	77	1.79
SLY9	$2.988^{{+}0.016}_{{-}0.008}$	0.86	$0.999^{{+}0.037}_{{-}0.019}$	0.23	—	91	1.68
BHF_BBB2	$3.021^{{+}0.020}_{{-}0.009}$	1.26	$2.032^{{+}0.053}_{{-}0.025}$	0.27	—	63	1.87
HQC18	$3.014^{{+}0.018}_{{-}0.010}$	1.03	$1.305^{{+}0.040}_{{-}0.020}$	0.22	—	81	1.72
KDE0V1	$3.005^{{+}0.018}_{{-}0.009}$	1.13	$1.837^{{+}0.048}_{{-}0.024}$	0.27	—	68	1.80
RS	$2.980^{{+}0.016}_{{-}0.008}$	0.80	$1.174^{{+}0.037}_{{-}0.019}$	0.24	—	92	1.61
SK272	$2.965^{{+}0.016}_{{-}0.008}$	0.75	$0.780^{{+}0.036}_{{-}0.019}$	0.23	No	100	1.43
SKI3	$2.967^{{+}0.015}_{{-}0.008}$	0.68	$0.744^{{+}0.034}_{{-}0.018}$	0.22	No	107	1.35
SKI5	$2.957^{{+}0.015}_{{-}0.008}$	0.63	$0.775^{{+}0.033}_{{-}0.019}$	0.23	No	109	1.32
SKMP	$2.990^{{+}0.016}_{{-}0.008}$	0.86	$1.186^{{+}0.037}_{{-}0.019}$	0.23	—	88	1.68
SLY	$3.005^{{+}0.017}_{{-}0.009}$	1.04	$1.425^{{+}0.042}_{{-}0.021}$	0.25	Possibly	76	1.82
SLY230A	$3.007^{{+}0.017}_{{-}0.009}$	0.97	$1.175^{{+}0.040}_{{-}0.020}$	0.23	—	83	1.76

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–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 multi-messenger constraints, see [119]). However, the theoretical modeling of the kilonova is a very complex and rapidly evolving field of research (e.g. [120–123] and the references therein).

The most relevant quantity regarding the remnant type is the baryonic mass $M_\mathrm{B}^\mathrm{R}$ of the system, in comparison to the maximum baryonic mass of non-rotating neutron stars, $\hat{M}_\mathrm{B}^\mathrm{TOV}$, and the maximum baryonic mass of neutron stars uniformly rotating at the mass-shedding limit, $\hat{M}_\mathrm{B}^\mathrm{ROT}$. Since the baryon number is conserved, the remnant mass is given by $M_\mathrm{B}^\mathrm{R} = M_\mathrm{B}^0 - M_\mathrm{B}^\mathrm{E}$, where $M_\mathrm{B}^\mathrm{E}$ denotes the baryonic mass expelled from the system dynamically and by winds, and $M_\mathrm{B}^0$ the total initial baryonic mass.

For discussing the expected remnant evolution, we introduce a useful dimensionless measure $\tilde{m}$ defined as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:def_mtilde} \tilde{m}(M_\mathrm{B}) = k_m \left(\frac{M_\mathrm{B}}{\hat{M}_\mathrm{B}^\mathrm{TOV}} - 1 \right), \qquad k_m = \frac{\hat{M}_\mathrm{B}^\mathrm{TOV}} {\hat{M}_\mathrm{B}^\mathrm{ROT} - \hat{M}_\mathrm{B}^\mathrm{TOV} }. \nonumber \end{align} \tag{ 4 }$$

The maximum masses $\hat{M}_\mathrm{B}^\mathrm{ROT}$ and $\hat{M}_\mathrm{B}^\mathrm{TOV}$, which we compute for each equation of state model are given in table A1. Based on $\newcommand{\e}{{\rm e}} \tilde{m}_\mathrm{R} \equiv \tilde{m}(M_\mathrm{B}^\mathrm{R})$ for a given equation of state model, we classify the remnant as hypermassive if $\tilde{m}_\mathrm{R}>1$, supramassive if $0 \leqslant \tilde{m}_\mathrm{R} \leqslant 1$, and stable if $\tilde{m}_\mathrm{R} < 0$.

The quantity $\tilde{m}_\mathrm{R}$ serves as a rough indication for the remnant lifetime $\tau_\mathrm{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 $\mathrm{ms}$. Within the narrow supramassive mass range, the lifetime decreases sharply. Computing $\tau_\mathrm{R}$ as a function of $\tilde{m}_\mathrm{R}$, however, is a difficult problem beyond the scope of this paper. Although we expect $\tilde{m}_\mathrm{R}$ to be the most important parameter, note that $\tau_\mathrm{R}$ might also depend on other factors such as mass ratio, spins, and equation of state.

In table 2, we provide the value $\tilde{m}_0$ corresponding to the median value of the total initial baryonic mass inferred for the narrow prior. Applying equation (4) to $M_\mathrm{B}^\mathrm{R}$ and $M_\mathrm{B}^0$, we obtain

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:remnant_mass} \tilde{m}_\mathrm{R} = \tilde{m}_0 - \left(\tilde{m}_0 + k_m \right) \frac{M_\mathrm{B}^\mathrm{E}}{M_\mathrm{B}^0}. \nonumber \end{align} \tag{ 5 }$$

In [116], the ejecta mass was estimated from the observed kilonova to around $0.078\, M_\odot$. If we assume more conservatively that $M_\mathrm{B}^\mathrm{E} < 0.1 \, M_\odot$, we find that $0 < \tilde{m}_0 - \tilde{m}_\mathrm{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 short-lived 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 $\tilde{m}_\mathrm{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_\odot$ due to the granularity of the simulated systems. The available thresholds are for equal-mass, non-spinning systems. To account for a possible impact of mass ratio and spin, we further increase the error by $0.05 \, M_\odot$. 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_\mathrm{B}^\mathrm{tot}$ 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_\mathrm{B}^0$ 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_\mathrm{B}^0$ are the result of a marginalization over mass ratio.

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) as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:max_tov_mass} \hat{M}_\mathrm{B}^\mathrm{TOV} = M_\mathrm{B}^0 \frac{1-M_\mathrm{B}^\mathrm{E} / M_\mathrm{B}^0}{1+ \tilde{m}_\mathrm{R} / k_m} \leqslant \frac{M_\mathrm{B}^0}{1+ \tilde{m}_\mathrm{R} / k_m}. \nonumber \end{align} \tag{ 6 }$$

From table 2, we find that $M_\mathrm{B}^0$ varies only slightly with equation of state. Therefore, we take the bound $M_\mathrm{B}^0 < 3.05 \,M_\odot$ 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 $\tilde{m}_\mathrm{R}$ the remnant lifetime is above the short gamma-ray burst delay. We already know that the lifetime is infinite unless $\tilde{m}_\mathrm{R}>0$. Future studies might provide tighter bounds, which will tighten the upper bound obtained for $\hat{M}_\mathrm{B}^\mathrm{TOV}$. For example, if we demand $\tilde{m}_\mathrm{R}>0.4$, then $\hat{M}_\mathrm{B}^\mathrm{TOV}<2.9 \, M_\odot$.

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 $\tilde{m}_\mathrm{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 $\tilde{m}_\mathrm{R} \geqslant 1$. From equation (6), we would obtain a corresponding limit $\hat{M}_\mathrm{B}^\mathrm{TOV}<2.67 \,M_\odot$ (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 $\hat{M}_\mathrm{G}^\mathrm{TOV} \leqslant 2.32 \,M_\odot$. 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 $\hat{M}_\mathrm{G}^\mathrm{TOV} \leqslant 2.28\pm 0.23\,M_\odot$ 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 $\hat{M}_\mathrm{G}^\mathrm{TOV} \leqslant 2.15$–$2.25 \, M_\odot$. It is based again on the detected limits for the total initial gravitational mass, and further assumes a total energy loss of $(0.15 \pm 0.03) \, M_\odot$, as well as a fixed difference $0.4\,M_\odot$ 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 $\hat{M}_\mathrm{G}^\mathrm{ROT} - \hat{M}_\mathrm{G}^\mathrm{TOV} = 0.34 \, M_\odot$ for the BHF_BBB2 equation of state). A recent analysis [125], which is not relying on assumptions about $\tilde{m}$, indicates that the remnant at the time of collapse might not rotate rapidly, and provides a limit $\hat{M}_\mathrm{G}^\mathrm{TOV} \leqslant 2.3 \, M_\odot$.

By using the detected gravitational-wave signal alone, without assumptions on the remnant fate, [40] inferred a limit $\hat{M}_\mathrm{G}^\mathrm{TOV} \leqslant 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.

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_\mathrm{B}$ and $\tilde{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 3. Implications derived assuming the MPA1 equation of state model and narrow prior, using different waveform approximants. The given quantities are explained in table 2.

Waveform	$M_\mathrm{B}$ $(M_\odot)$	$E_c / c^2$ $(10^{15}\, \mathrm{g}~\mathrm{cm}^{-3})$	$\tilde{m}_0$	$\Delta_{0.1}$	$I_\mathrm{max}$ $(M_\odot^3)$	$F_\mathrm{max}$ $(\mathrm{kHz})$
TaylorF2(7.5 PN)	$3.004^{{+}0.016}_{{-}0.009}$	0.74	$-0.016^{{+}0.029}_{{-}0.016}$	0.18	114	1.34
TaylorF2(7 PN)	$3.004^{{+}0.016}_{{-}0.009}$	0.74	$-0.017^{{+}0.029}_{{-}0.015}$	0.18	114	1.34
TaylorF2(6 PN)	$3.004^{{+}0.016}_{{-}0.009}$	0.74	$-0.016^{{+}0.028}_{{-}0.016}$	0.18	114	1.34
EOB+NRT	$3.004^{{+}0.016}_{{-}0.008}$	0.74	$-0.015^{{+}0.027}_{{-}0.015}$	0.18	114	1.34
Phenom+NRT	$3.005^{{+}0.016}_{{-}0.009}$	0.74	$-0.015^{{+}0.029}_{{-}0.015}$	0.18	114	1.34

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.

Table 4. Like table 2, but showing results for the broad prior, again with the TaylorF2 (7.5 PN) waveform model.

Equation of state	$M_\mathrm{B}^0$ $(M_\odot)$	$E_c / c^2$ $(10^{15}\, \mathrm{g}~\mathrm{cm}^{-3})$	$\tilde{m}_0$	$\Delta_{0.1}$	Prompt collapse	$I_\mathrm{max}$ $(M_\odot^3)$	$F_\mathrm{max}$ $(\mathrm{kHz})$
APR4_EPP	$3.084^{{+}0.316}_{{-}0.070}$	1.48	$1.019^{{+}0.683}_{{-}0.150}$	0.22	Possibly	84	1.89
MS1_PP	$2.960^{{+}0.121}_{{-}0.026}$	0.53	$-0.574^{{+}0.202}_{{-}0.044}$	0.17	No	158	1.03
MS1B_PP	$2.973^{{+}0.144}_{{-}0.030}$	0.55	$-0.542^{{+}0.231}_{{-}0.048}$	0.16	No	159	1.05
KDE0V	$3.058^{{+}0.219}_{{-}0.053}$	1.80	$1.978^{{+}0.582}_{{-}0.139}$	0.27	—	66	1.86
SKI6	$3.028^{{+}0.225}_{{-}0.041}$	1.08	$0.922^{{+}0.489}_{{-}0.087}$	0.22	—	97	1.66
H4	$3.007^{{+}0.178}_{{-}0.035}$	0.94	$1.450^{{+}0.388}_{{-}0.077}$	0.22	Possibly	104	1.36
MPA1	$3.041^{{+}0.255}_{{-}0.043}$	0.90	$0.050^{{+}0.456}_{{-}0.077}$	0.18	No	121	1.41
SK255	$3.003^{{+}0.203}_{{-}0.041}$	1.10	$1.207^{{+}0.497}_{{-}0.102}$	0.24	—	91	1.61
SKI2	$3.001^{{+}0.176}_{{-}0.038}$	0.97	$1.100^{{+}0.408}_{{-}0.089}$	0.23	—	100	1.53
SKI4	$3.028^{{+}0.228}_{{-}0.037}$	1.11	$0.978^{{+}0.494}_{{-}0.081}$	0.22	—	95	1.67
SKOP	$3.019^{{+}0.190}_{{-}0.032}$	1.45	$1.874^{{+}0.500}_{{-}0.085}$	0.26	—	73	1.70
SLY2	$3.043^{{+}0.239}_{{-}0.045}$	1.48	$1.494^{{+}0.585}_{{-}0.110}$	0.25	—	77	1.79
SLY9	$3.018^{{+}0.234}_{{-}0.037}$	1.15	$1.068^{{+}0.537}_{{-}0.083}$	0.23	—	91	1.68
BHF_BBB2	$3.077^{{+}0.201}_{{-}0.063}$	2.02	$2.181^{{+}0.545}_{{-}0.168}$	0.27	—	63	1.87
HQC18	$3.064^{{+}0.253}_{{-}0.057}$	1.42	$1.418^{{+}0.561}_{{-}0.129}$	0.22	—	81	1.72
KDE0V1	$3.049^{{+}0.216}_{{-}0.051}$	1.73	$1.953^{{+}0.580}_{{-}0.135}$	0.27	—	68	1.80
RS	$3.016^{{+}0.216}_{{-}0.042}$	1.14	$1.259^{{+}0.507}_{{-}0.099}$	0.24	—	92	1.61
SK272	$3.003^{{+}0.209}_{{-}0.044}$	1.00	$0.868^{{+}0.483}_{{-}0.102}$	0.23	—	100	1.61
SKI3	$3.000^{{+}0.190}_{{-}0.039}$	0.91	$0.818^{{+}0.421}_{{-}0.087}$	0.22	—	107	1.53
SKI5	$2.975^{{+}0.116}_{{-}0.024}$	0.78	$0.816^{{+}0.262}_{{-}0.056}$	0.23	—	109	1.47
SKMP	$3.019^{{+}0.216}_{{-}0.035}$	1.19	$1.254^{{+}0.504}_{{-}0.083}$	0.23	—	88	1.68
SLY	$3.042^{{+}0.238}_{{-}0.044}$	1.50	$1.517^{{+}0.590}_{{-}0.109}$	0.25	Possibly	76	1.82
SLY230A	$3.044^{{+}0.251}_{{-}0.044}$	1.36	$1.261^{{+}0.581}_{{-}0.101}$	0.23	—	83	1.76

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 $\tilde{m}_\mathrm{R}\geqslant 0$, we obtain $\hat{M}_\mathrm{B}^\mathrm{TOV} < 3.4 M_\odot$. Assuming $\tilde{m}_\mathrm{R}\geqslant1$, we would obtain $\hat{M}_\mathrm{B}^\mathrm{TOV} < 2.98 M_\odot$.

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\times 10^{15}~\mathrm{g}~\mathrm{cm}^{-3}$ (compare table 4).

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.