Inferring prompt black-hole formation in neutron star mergers from gravitational-wave data

TL;DR

This paper develops two methods to infer the likelihood of prompt black hole formation in neutron star mergers from gravitational-wave data, combining data analysis with numerical relativity insights.

Contribution

It introduces two novel approaches to estimate prompt black hole formation probability using inspiral gravitational-wave signals and equation-of-state models.

Findings

Both methods predict a 50-70% probability of prompt black hole formation for GW170817.

Probability estimates decrease below 10% when maximum mass constraints are applied.

The methods are consistent and provide insights into merger outcomes based on gravitational-wave data.

Abstract

The gravitational-wave GW170817 is associated to the inspiral phase of a binary neutron star coalescence event. The LIGO-Virgo detectors sensitivity at high frequencies was not sufficient to detect the signal corresponding to the merger and post-merger phases. Hence, the question whether the merger outcome was a prompt black hole formation or not must be answered using either the pre-merger gravitational wave signal or electromagnetic counterparts. In this work we present two methods to infer the probability of prompt black hole formation, using the analysis of the inspiral gravitational-wave signal. Both methods combine the posterior distribution from the gravitational-wave data analysis with numerical relativity results. One method relies on the use of phenomenological models for the equation of state and on the estimate of the collapse threshold mass. The other is based on the estimate of the tidal polarizability parameter $\tilde{\Lambda}$ that is correlated in an equation-of-state agnostic way with the prompt BH formation. We analyze GW170817 data and find that the two methods consistently predict a probability of ~ 50-70% for prompt black-hole formation, which however may significantly decrease below 10% if the maximum mass constraint from PSR J0348+0432 or PSR J0740+6620 is imposed.

Tables

Table 1. Table 1: Summary of injections and TaylorF2 recovery. Collapse time t BH subscript 𝑡 BH t_{\text{BH}} is reported from merger time, defined at the peak of the amplitude. We indicate with HMNS (MNS) remnants that are short (long) lived, i.e. that (do not) collapse to BH within the simulated time. All the simulations are performed at standard resolution of Radice  et al. ( 2018b ) .

EOS	$M_{\text{TOV}}^{\text{max}}$	$C_{\text{max}}$	$M_{\text{thr}}$	$M_A$	$M_B$	${ℳ}_c$	${\Lambda }_2^A$	${\Lambda }_2^B$	$\tilde{\Lambda }$	$t_{\text{BH}}$	Remnant at	Ref.	$P_{\mathrm{PC}}^{M_{\text{thr}}}$	$P_{\mathrm{PC}}^{M_{\text{thr}},M_{\text{TOV}}^{\text{max}}}$	$P_{\mathrm{PC}}^{{\tilde{\Lambda }}_{\mathrm{thr}}}$
	$\left[M_{\odot }\right]$		$\left[M_{\odot }\right]$	$\left[M_{\odot }\right]$	$\left[M_{\odot }\right]$	$\left[M_{\odot }\right]$				[ms]	$t\sim 3$ ms		%	%	%
2B	${1.78}^{\ast }$	0.3120	2.43${}_{-0.24}{}^{}{}_{}{}^{+0.24}$	1.35	1.35	1.17	127	127	127	0.49	BH	Bernuzzi et al. (2015)	$99.5$	${48.3}^{\ast }$	$100$
SLy	2.06	0.3066	2.87${}_{-0.30}{}^{}{}_{}{}^{+0.26}$	1.50	1.50	1.30	191	191	191	0.99	BH	Dietrich et al. (2018)	$94.5$	$82.9$	$100$
LS220	2.04	0.2841	2.95${}_{-0.24}{}^{}{}_{}{}^{+0.29}$	1.60	1.60	1.39	202	202	202	0.63	BH	Radice et al. (2018c, a, b)	$90.0$	$84.4$	$99.9$
SFHo	2.06	0.2952	2.95${}_{-0.29}{}^{}{}_{}{}^{+0.25}$	1.46	1.46	1.27	252	252	252	0.70	BH	Radice et al. (2018c, a, b)	$72.7$	$40.4$	$97.8$
BHB$\Lambda \varphi$	2.11	0.2677	3.10${}_{-0.18}{}^{}{}_{}{}^{+0.35}$	1.60	1.60	1.39	306	306	306	0.99	BH	Radice et al. (2017, 2018c, 2018a, 2018b)	$36.0$	$19.6$	$72.9$
DD2	2.42	0.3007	3.35${}_{-0.28}{}^{}{}_{}{}^{+0.37}$	1.59	1.59	1.38	332	332	332	$\sim 3$	BH	Radice et al. (2017)	$32.1$	$18.5$	$66.2$
SFHo	2.06	0.2952	2.95${}_{-0.29}{}^{}{}_{}{}^{+0.25}$	1.40	1.40	1.22	334	334	334	1.07	BH	Radice et al. (2018c, a)	$41.7$	$4.9$	$60.6$
ALF2	1.99	0.2602	2.87${}_{-0.06}{}^{}{}_{}{}^{+0.43}$	1.50	1.50	1.30	382	382	382	0.64	BH	Dietrich et al. (2018)	$20.9$	$3.6$	$29.2$
SLy-SOR	2.06	0.3066	2.87${}_{-0.30}{}^{}{}_{}{}^{+0.26}$	1.34	1.34	1.17	401	401	401	$\sim 14$	HMNS	This work	$25.0$	$0.4$	$21.8$
SLy-SOR	2.06	0.3066	2.87${}_{-0.30}{}^{}{}_{}{}^{+0.26}$	1.43	1.26	1.17	264	592	401	$\sim 13$	HMNS	This work	$23.7$	$0.1$	$21.6$
SFHo	2.06	0.2952	2.95${}_{-0.29}{}^{}{}_{}{}^{+0.25}$	1.44	1.27	1.18	274	606	412	$\sim 12$	HMNS	This work	$21.6$	$2.0$	$17$
SFHo	2.06	0.2952	2.95${}_{-0.29}{}^{}{}_{}{}^{+0.25}$	1.35	1.35	1.18	413	413	413	$\sim 4$	HMNS	Radice et al. (2018c, a, b)	$20.5$	$0.2$	$15.6$
LS220	2.04	0.2841	2.95${}_{-0.24}{}^{}{}_{}{}^{+0.29}$	1.44	1.25	1.17	432	1136	713	$\sim 33$	HMNS	This work	$2.4$	$0.0$	$0$
LS220	2.04	0.2841	2.95${}_{-0.24}{}^{}{}_{}{}^{+0.29}$	1.34	1.34	1.17	715	715	715	$\sim 16$	HMNS	This work	$0.7$	$0.0$	$0$
DD2	2.42	0.3007	3.35${}_{-0.28}{}^{}{}_{}{}^{+0.37}$	1.36	1.36	1.18	840	840	840	$\sim 21$	MNS	Perego et al. (2019)	$0.1$	$0.0$	$0$
DD2	2.42	0.3007	3.35${}_{-0.28}{}^{}{}_{}{}^{+0.37}$	1.24	1.24	1.08	1366	1366	1366	$>20$	MNS	Radice et al. (2018a)	$0.0$	$0.0$	$0$
BHB$\Lambda \varphi$	2.11	0.2677	3.10${}_{-0.18}{}^{}{}_{}{}^{+0.35}$	1.24	1.24	1.08	1367	1367	1367	$>20$	MNS	Radice et al. (2017)	$0.0$	$0.0$	$0$

Table 2. Table 2: Summary of GW170817 results derived with the threshold-mass and the threshold- Λ ~ ~ Λ \tilde{\Lambda} methods: probability of prompt collapse for different sets of analyses published by the LVC.

Method	Inferred parameters	Approximant	Ref.	$P_{\text{PC}}$
$M_{\text{thr}}$	$\overrightarrow{\gamma }$	IMRPhenomPv2NRtidal	Abbott et al. (2018b)	0.59
$M_{\text{thr}}$	$\overrightarrow{\gamma }|M_{\text{TOV}}^{\text{max}}\ge 1.97M_{\odot }$	IMRPhenomPv2NRtidal	Abbott et al. (2018b)	0.09
${\tilde{\Lambda }}_{\text{thr}}$	$\overrightarrow{\gamma }$	IMRPhenomPv2NRtidal	Abbott et al. (2018b)	0.69
${\tilde{\Lambda }}_{\text{thr}}$	$\overrightarrow{\gamma }|M_{\text{TOV}}^{\text{max}}\ge 1.97M_{\odot }$	IMRPhenomPv2NRtidal	Abbott et al. (2018b)	0.44
${\tilde{\Lambda }}_{\text{thr}}$	$\tilde{\Lambda }$	TaylorF2	Abbott et al. (2019)	0.54
${\tilde{\Lambda }}_{\text{thr}}$	$\tilde{\Lambda }$	IMRPhenomPv2NRtidal	Abbott et al. (2019)	0.58
${\tilde{\Lambda }}_{\text{thr}}$	$\tilde{\Lambda }$	IMRPhenomDNRtidal	Abbott et al. (2019)	0.59
${\tilde{\Lambda }}_{\text{thr}}$	$\tilde{\Lambda }$	SEOBNRT	Abbott et al. (2019)	0.60
${\tilde{\Lambda }}_{\text{thr}}$	${\Lambda }_s$	IMRPhenomPv2NRtidal	Abbott et al. (2018b)	0.74

Table 3. Table 3: Numerical relativity data used for determinig the prompt collapse mass threshold coefficient k thr subscript 𝑘 thr k_{\text{thr}} .

EOS	$C_{\text{max}}$	$k_{\text{thr}}$	$\delta k_{\text{thr}}$	$M_{\text{TOV}}^{\text{max}}$	$M_{\text{thr}}$	Ref
				$\left[M_{\odot }\right]$	$\left[M_{\odot }\right]$
APR4	0.329	1.243	0.023	2.21	2.80	Hotokezaka et al. (2011)
SLy	0.307	1.342	0.024	2.06	2.76	Hotokezaka et al. (2011)
H3	0.224	1.566	0.056	1.79	2.90	Hotokezaka et al. (2011)
H4	0.258	1.452	0.025	2.03	2.95	Hotokezaka et al. (2011)
ALF2	0.260	1.414	0.024	1.99	2.81	Hotokezaka et al. (2011)
NL3	0.307	1.380	0.018	2.79	3.85	Bauswein et al. (2013)
GS1	0.306	1.400	0.018	2.75	3.85	Bauswein et al. (2013)
LS375	0.325	1.347	0.018	2.71	3.65	Bauswein et al. (2013)
DD2	0.300	1.384	0.021	2.42	3.35	Bauswein et al. (2013)
Shen	0.250	1.554	0.023	2.22	3.45	Bauswein et al. (2013)
TM1	0.260	1.561	0.023	2.21	3.45	Bauswein et al. (2013)
SFHx	0.292	1.432	0.023	2.13	3.05	Bauswein et al. (2013)
GS2	0.262	1.555	0.024	2.09	3.25	Bauswein et al. (2013)
SFHo	0.294	1.432	0.024	2.06	2.95	Bauswein et al. (2013)
LS220	0.284	1.495	0.025	2.04	3.05	Bauswein et al. (2013)
TMA	0.247	1.609	0.025	2.02	3.25	Bauswein et al. (2013)
IUF	0.255	1.564	0.026	1.95	3.05	Bauswein et al. (2013)
LS220	0.284	1.445	0.024	2.04	2.95	Zappa et al. (2018)
BHB$\Lambda \varphi$	0.268	1.469	0.047	2.11	3.10	Zappa et al. (2018)
ALF2	0.260	1.444	0.063	1.99	2.86	Zappa et al. (2018)
H4	0.258	1.529	0.049	2.03	3.10	Zappa et al. (2018)
SLy	0.307	1.395	0.061	2.06	2.86	Zappa et al. (2018)
BHB$\Lambda \varphi$	0.268	1.503	0.024	2.11	3.17	Köppel et al. (2019)
DD2	0.300	1.364	0.021	2.42	3.30	Köppel et al. (2019)
SFHo	0.294	1.391	0.024	2.06	2.87	Köppel et al. (2019)
TM1	0.260	1.520	0.023	2.21	3.36	Köppel et al. (2019)

Equations

$$M_{\text{thr}}=k_{\text{thr}}{M}_{\text{TOV}}^{\text{max}}.$$

$$1.3\lesssim k_{\text{thr}}\lesssim 1.7.$$

$$k_{\text{thr}}(C_{\text{max}})=-(3.29\pm 0.23)C_{\text{max}}+(2.392\pm 0.064).$$

$${\kappa }_2^{\text{T}}$$

$${\Lambda }_2^{\text{A}}=\frac{2}{3}{k}_2^{\text{A}}{\left(\frac{c^2}{G}\frac{R_{\text{A}}}{M_{\text{A}}}\right)}^5,$$

$$\tilde{\Lambda }$$

$$338\lesssim {\tilde{\Lambda }}_{\text{thr}}\lesssim 386.$$

$$\Gamma =\exp \left[\sum _{k=0}^3{\gamma }_k\log (p/p_0{)}^k\right].$$

$$P_{\mathrm{PC}}=P(M>M_{\text{thr}}(\vec{\gamma })|d),$$

$$P(\text{prompt collapse}|\tilde{\Lambda })=\frac{1}{1+e^{\frac{\tilde{\Lambda }-{\tilde{\Lambda }}_0}{\beta }}},$$

$$P_{\mathrm{PC}}=\int \mathrm{d}\tilde{\Lambda }P(\text{prompt collapse}|\tilde{\Lambda })p(\tilde{\Lambda }|d).$$

$$w(\vec{\gamma })=p(M_{\text{TOV}}^{\text{max}}(\vec{\gamma })>\overline{M_{\max }}).$$

Full text

I Introduction

The gravitational-wave (GW) signal GW170817, detected by the LIGO-Virgo detector network Aasi et al. (2015); Acernese et al. (2015), is a chirp transient compatible with the emission from a binary neutron star system coalescence in the late-inspiral phase Abbott et al. (2017a, 2019, 2018a). The signal has significant signal-to-noise ratio (SNR) in the range 30 to 600 Hz, roughly corresponding to the last 100 to 30 orbits to merger for an equal-mass binary with total mass $M\sim 2.7\>{\rm M_{\odot}}$. The data analysis of GW170817 provided us with an estimate of the dominant tidal polarizability parameter that, in turn, constrains the NS cold equation of state Damour et al. (2012); Favata (2014); Del Pozzo et al. (2013); Lackey and Wade (2015); Abbott et al. (2018b). The LIGO-Virgo detectors’ sensitivity was not sufficient to detect a signal from the merger phase and the remnant, which lie in the kHz range Abbott et al. (2018c). An outstanding question is thus whether the coalescence resulted in the formation of a black-hole (BH) or in a NS remnant.

A first answer was given by the interpretation of the electromagnetic counterparts observed with delays of seconds to days with respect to the GW and composed by a GRB Abbott et al. (2017b); Troja et al. (2017); Abbott et al. (2017c) and a kilonova Coulter et al. (2017); Chornock et al. (2017); Nicholl et al. (2017); Cowperthwaite et al. (2017); Tanvir et al. (2017); Tanaka et al. (2017). Energetics and timing of the latter exclude both a prompt BH formation and a long-lived remnant, e.g. Margalit and Metzger (2017); Shibata et al. (2017). Most likely, the merger dynamics produced a hypermassive NS that collapsed on timescales of ${\sim}0.01$ to ${\sim}2$ seconds. Such a conclusion is informed and supported by numerical relativity (NR) results that established the formation of hypermassive NS remnants for canonical NS masses and equations of state supporting $M_{\text{TOV}}^{\text{max}}\gtrsim 2{\rm M_{\odot}}$ Shibata and Uryu (2000); Shibata et al. (2003); Anderson et al. (2008); Baiotti et al. (2008); Thierfelder et al. (2011); Hotokezaka et al. (2011); Bauswein et al. (2013); Shibata et al. (2017); Radice et al. (2018a).

In this work, we explore a different approach to inferring the merger remnant. Instead of considering the EM counterparts, we consider the pre-merger GW and infer binary parameters using the late-inspiral solely. The posterior distributions of these parameters are then combined with information from NR simulations. Our methods allow us to quantify the probability that a BH was promptly formed.

This paper is structured as follows: Sec. II outlines the input from NR data in our inference methods, based on which we classify the outcome of a BNS merger; the two methods are introduced in Sec. III and are validated by analyzing a set of simulated GW detections in Sec. IV; we perform the analysis on GW170817 data and present our results in Sec. V, while some concluding remarks are given in Sec. VI. We use geometric units $G=c=1$ unless stated differently.

II Prompt Collapse Threshold

II.1 Mass threshold estimate

Numerical-relativity simulations indicate that a NS binary merger will be followed by a prompt collapse to a BH, if the total gravitational mass $M$ of the binary exceeds a threshold mass, given by Hotokezaka et al. (2011); Bauswein et al. (2013)

[TABLE]

In the expression above, $k_{\text{thr}}$ depends, in general, on the EOS, mass ratio, and spin, while $M_{\text{TOV}}^{\text{max}}$ is the gravitational mass of the heaviest stable nonrotating NS, which also depends on the EOS. Empirically, the prompt collapse threshold is calculated from the simulations by considering remnants that collapse within 2 ms from the waveform peak amplitude (conventionally, the “merger time”). Examples of merger waveforms for a prompt collapse and a NS remnant are shown in Fig. 1. In the prompt collapse case an apparent horizon forms during the simulation at a time close to the retarded merger time; the waveform frequency at those times corresponds to the quasi-normal mode of the black hole.

For a sample of hadronic EOS and equal-mass nonspinning binaries, the threshold coefficient in Eq. (1) is found in the range Hotokezaka et al. (2011); Bauswein et al. (2013); Köppel et al. (2019)

[TABLE]

Considering a sample of equal-mass, nonspinning binaries and 12 hadronic EOS, Ref. Bauswein et al. (2013) showed that $k_{\text{trh}}$ has an approximately EOS-independent linear behavior in the maximum compactness $C_{\text{max}}$ of nonrotating equilibrium NS solution. Note that by inverting Eq. (1) and assuming that the merger did not promptly form a BH, one may obtain a bound on the maximum stable NS mass Shibata et al. (2017); Bauswein et al. (2017); Ruiz et al. (2018).

We have repeated the analysis on the threshold mass with the data of CoRe collaboration Zappa et al. (2018); Dietrich et al. (2018) by including 10 new simulations with 5 EOS and different masses and spins. We have compared and combined our new results with the ones reported in Hotokezaka et al. (2011); Bauswein et al. (2013); Köppel et al. (2019). Our final sample includes 18 different hadronic EOS and for 8 of them results from more than one analysis are available. Using the results reported by Bauswein et al. (2013) and Köppel et al. (2019), and by adding the data of CoRe collaboration, we find a linear fit with updated coefficients that reads

[TABLE]

The data that were used for this fit are given in Appendix A, along with further details.

II.2 Tidal parameter threshold estimate

Alternatively, the prompt collapse threshold can be characterized in an EOS-independent way in terms of the tidal polarizability parameter

[TABLE]

where the tidal polarizability coefficient of star $A$ is

[TABLE]

and $k^{\text{A}}_{2}$ is the quadrupolar gravito-electric Love number Hinderer (2008); Damour and Nagar (2009); Binnington and Poisson (2009). Above, $(R_{\text{A}},M_{\text{A}})$ are the NS areal radius and mass and $X_{\text{A}}=M_{\text{A}}/M$. The $\Lambda_{2}$ parameter is strongly dependent on the NS internal structure; thus, its measurement provides a constraint on the NS EOS 111Black holes are not deformed in this way; black hole static perturbations lead to $k_{2}=0$ Damour and Lecian (2009); Damour and Nagar (2010); Binnington and Poisson (2009); Gürlebeck (2015).. The binary’s post-Newtonian tidal dynamics and waveform are parametrized at leading-order by $\kappa^{\text{T}}_{2}$ Damour et al. (2012); Bernuzzi et al. (2014). A tidal polarizability parameter commonly used in GW analysis (and equivalent to $\kappa^{\text{T}}_{2}$ for equal-mass binaries) is

[TABLE]

By analyzing the NR data of the CoRe collaboration, Ref. Zappa et al. (2018) found that all the reported prompt collapse mergers are captured by the condition $\kappa^{\text{T}}_{2}<73$, with a variability of $\delta\kappa^{\text{T}}_{2}<40$, depending on the EOS. Inspection of the same NR data provided also a range for the value of $\tilde{\Lambda}$ at the prompt collapse threshold:

[TABLE]

III Method

Based on the universal behavior discussed in Sec. II, we present two different ways of inferring whether a BNS merger is followed by a prompt collapse to a BH using solely GW data (with the exception of the sky location which we may fix to the one obtained by EM observations, when an EM counterpart is available). We test the validity of our methods against a set of high-resolution numerical simulations of BNS mergers with different masses and EOS.

For our Bayesian data analysis on the GW signal, we use a Markov-chain Monte Carlo (MCMC) algorithm as implemented in the LALInference software package Veitch et al. (2015), with a set-up similar to the one employed in the latest LVC analysis of GW170817 Abbott et al. (2019, 2018b).

III.1 Threshold Mass

For this method we make use of the mass threshold estimate of Sec. II.1, whereby the total mass $M$ of the progenitor NS binary being larger or smaller than $M_{\text{thr}}$ determines whether the product of the merger will promptly collapse to a BH or not. The threshold mass $M_{\text{thr}}$ depends on the EOS via Eq. (1) and (3). We perform a full Bayesian analysis on the data, that returns posterior distributions for the binary parameters, including the EOS. The barotropic EOS for the cold dense NS matter is sampled through a 4-dimensional family of pressure-density functions $P(\rho)$, parametrized by $(\gamma_{0},\gamma_{1},\gamma_{2},\gamma_{3})$ in the spectral decomposition Lindblom (2010); Carney et al. (2018); Abbott et al. (2018b), a smooth alternative to piecewise polytropic models Read et al. (2009); Lackey and Wade (2015), where the adiabatic index $\Gamma=\rho\frac{{\rm d}\ln P}{{\rm d}\rho}$ is given by

[TABLE]

with $p_{0}$ some reference pressure. For each sampled point in the parameter space, we solve the TOV equations to calculate not only the tidal polarizability parameters $\Lambda_{2}^{A}$ which are used to model the tidal effects in the waveform, but also the values for $M_{\text{TOV}}^{\text{max}}$, $C_{\rm max}$ and $k_{\text{thr}}$. We can thus translate the joint posterior PDF on masses and EOS parameters $(m_{1},m_{2},\gamma_{0},\gamma_{1},\gamma_{2},\gamma_{3})$ into a joint posterior PDF on the $(M,M_{\text{thr}})$ plane. The fraction of the posterior distribution that lies above the diagonal is equal to the posterior probability of prompt collapse

[TABLE]

where $\vec{\gamma}=(\gamma_{0},\gamma_{1},\gamma_{2},\gamma_{3})$ and $d$ denotes our data.

As an additional step, one may choose to impose further implicit constraints on the parameter space, such as requiring that the EOS support NS masses larger than a given value. For instance, the observation of the binary pulsar PSR J0348+0432 Antoniadis et al. (2013) and the narrow measurement of the pulsar’s mass, or even the more recent measurement of an even heavier (but with larger uncertainty) pulsar mass in J0740+6620 Cromartie et al. (2019). In one of the analyses of Abbott et al. (2018b) the conservative 1-$\sigma$ bound for PSR J0348+0432 at $1.97{\rm M_{\odot}}$ was considered as a hard constraint. Here we take a different approach and marginalize over the mass measurement uncertainties into our analysis by treating that measurement as a random variable sampled from a normal distribution that is adapted to the mean and standard deviation of the measurement, and weighing the posterior samples accordingly.

III.2 Threshold Tidal Parameter

For the second method we again employ a Bayesian analysis of the GW data, this time focusing on the posterior distribution of the tidal deformability parameter $\tilde{\Lambda}$ given by Eq. (6). The set-up of our Bayesian analysis follows that of Abbott et al. (2019). We then make use of the corresponding criterion of Eq. (7) in order to estimate the probability of prompt collapse. Note that the criterion defines a transition region between the studied cases where the merger product undergoes prompt collapse and the ones where it does not. The outcomes of NR simulations within this transition region are not perfectly ordered. We treat this classification problem by assigning a probability distribution to the uncertainty of the threshold value $\tilde{\Lambda}_{\text{thr}}$ instead of choosing a hard threshold or, equivalently, by defining a sigmoid-type conditional probability of prompt collapse for a given value of $\tilde{\Lambda}$ as

[TABLE]

which tends to $1$ ([math]) for small (large) values of $\tilde{\Lambda}$. The values for the sigmoid parameters, i.e. the central value and the width, are chosen based on the available set of NR simulations in this region to be $\tilde{\Lambda}_{0}\approx 362$ and $\beta\approx 13.7$ respectively. Then, once the posterior PDF $p(\tilde{\Lambda}|d)$ is calculated, the probability of prompt collapse is simply computed by integrating the posteriors from the minimum value up to the threshold value using the sigmoid of Eq. (10) as a kernel

[TABLE]

Note that this method does not rely on any assumption about the EOS, but only on the phenomenological parameter $\tilde{\Lambda}$ which is directly measured from the data. In the present analysis we assume that Eq. (7) holds independently on $q$ and spins. That hypothesis is justified by inspection of the CoRe data that span $q\in[1,2]$ and dimensionless spins up to $\sim 0.1$.

IV Injection studies

We validate our methods using injections of known inspiral-merger waveforms corresponding to binaries simulated in NR. We demonstrate that both methods are effective in estimating the collapse treshold and discuss/quantify their systematics.

IV.1 Setup

We consider NR merger simulations of irrotational binaries with different chirp masses specifically performed for this work together with data previously presented in Bernuzzi et al. (2015); Radice et al. (2017, 2018c); Dietrich et al. (2018); Radice et al. (2018a, b); Perego et al. (2019). The new simulations are performed with the WhiskyTHC code Radice and Rezzolla (2012); Radice et al. (2014a, b) at multiple grid resolutions, using the same setup described in Radice et al. (2018c, a, b).

The main properties of the simulated binaries, the outcome of the merger and the summary data from the injection are summarized in Tab. 1. We simulate with five microphysical EOS: the BHB$\Lambda\phi$ EOS Banik et al. (2014), the DD2 EOS Typel et al. (2010); Hempel and Schaffner-Bielich (2010), the LS220 EOS Lattimer and Swesty (1991), the SFHo EOS Steiner et al. (2013), the SLy-SOR EOS Schneider et al. (2017); and tthree piecewise polytropic: the ALF2, 2B and the SLy EOS Read et al. (2009). The microphysical EOS predict NS maximum masses and radii within the range allowed by current astrophysical constraints. The 2B EOS is representative of soft EOS that do not support the largest NS masses observed so far Antoniadis et al. (2013); Cromartie et al. (2019). All the simulations with chirp mass $\mathcal{M}_{c}\sim 1.18\>{\rm M_{\odot}}$ and tidal deformability compatible with GW170817 except 2B, predict a short-lived NS remnant collapsing to BH within ${\sim}15$ ms. The simulation DD2 1.59+1.59 with $\tilde{\Lambda}=332$ is below the $\tilde{\Lambda}$ threshold for prompt BH formation but forms a short lived NS with lifetime ${\sim}3$ ms Radice et al. (2017). While this is possibly related to numerical uncertainties, the binary provides an interest borderline case.

Throughout this work we use the terms hypermassive NS (HMNS) and massive NS (MNS) with a slightly different meaning than what is commonly used in the literature 222A HMNS is defined as a differentially rotating NS at equilbrium with mass above the uniformly rotating limit Baumgarte et al. (2000). A supramassive NS is a rotating NS at equilibrium with rest mass exceeding the nonrotating limit $M_{\text{TOV}}^{\text{max}}$ Cook et al. (1994). A remnant with mass below $M_{\text{TOV}}^{\text{max}}$ is usually indicated as MNS. We indicate with HMNS (MNS) merger remnants that are short (long) lived, i.e. that (do not) collapse to BH within the simulated time. The reason for this choice is that merger remnants are not cold equilibrium NS configurations, and their secular evolution is far from being understood (see e.g. discussion in Radice et al. (2018a).)

The simulations provide us with dynamics and waveform starting from GW frequencies $\sim\,500-900$ Hz, depending on the binary mass and simulation length. Hence, the NR waveform alone are not sufficient to perform injection of BNS signals. Waveforms spanning from an initial GW frequency of 30 Hz to merger and corresponding to the binaries of Tab. 1 are constructed with the TEOBResumS waveform model Bernuzzi et al. (2015); Nagar et al. (2018). Specifically, we use the nonspinning tidal model of Akcay et al. (2019) with gravitational-self-force resummed gravitoelectric and post-Newtonian gravitomagnetic terms (Model GSF23*(+)PN(-)* with $p=4$ of Tab. I in Akcay et al. (2019)). Waveforms are generated using the post-adiabatic inspiral speed-up developed in Nagar and Rettegno (2019) 333The public available TEOBResumS code can be found at teo . .

For our Bayesian data analysis on the simulated GW signals, we use a Markov-chain Monte Carlo (MCMC) algorithm as implemented in the LALInference software package Veitch et al. (2015); LIGO Scientific Collaboration (2018), with a set-up similar to the one employed in the latest LVC analysis of GW170817 Abbott et al. (2019, 2018b). The simulated signals are coherently projected and analyzed as the output strain of LIGO Handford (H1), LIGO Livingston (L1) and Virgo (V1) at design sensitivity. The intrinsic parameters of the nonspinning BNS sources are given in Tab. 1, while the location and orientation parameters are compatible with GW170817. In order to isolate possible systematics from statistical uncertainties due to noise, we perform our tests on the zero-noise realization.

We perform our analyses using two different waveform models, namely TaylorF2 and IMRPhenomPv2NRtidal, both restricted to aligned dimensionless spins ranging within $[-0.05,0.05]$ (the low-spin prior of Abbott et al. (2019)). The tidal effects are modeled in the case of TaylorF2 up to 2.5PN beyond leading order Damour et al. (2012) and in the case of IMRPhenomPv2 using the NRTidal approach of Dietrich et al. (2017).

In the threshold-mass method, we are able to additionally impose an observational constraint on the EOS prior, coming from the heaviest observed NS. This can be either a hard constraint at a chosen mass value (e.g. $1.97$ ${\rm M_{\odot}}$ as in Abbott et al. (2018b)), or a probabilistic constraint that takes into account the measured posterior PDF. In the latter case we will use the median and 1-$\sigma$ error of the mass measurement of PSR J0348+0432 Antoniadis et al. (2013) to reconstruct a Gaussian PDF for the heaviest observed NS mass $p(\overline{M_{\rm max}})=\mathcal{N}(2.01{\rm M_{\odot}},0.04{\rm M_{\odot}})$ and assign a weight $w$ on each sampled EOS appropriately, according to its maximum NS mass,

[TABLE]

A comparison between results derived with and without such a constraint is demonstrated in Appendix C.

IV.2 Results

We find that the IMRPhenomPv2NRtidal waveform systematically underestimates the inference of the injected $\tilde{\Lambda}$. This result was anticipated by the high SNR injections of Dudi et al. (2018); Messina et al. (2019), but could not be studied systematically due to the limited number of injections performed there. A similar bias is present in the EOS inference runs with IMRPhenomPv2NRtidal, but the mass threshold method to determine the prompt collapse is less affected by systematics on tidal parameters than the $\tilde{\Lambda}$ threshold method. In the following, we discuss only the results obtained with TaylorF2. The effect of waveform systematics on the results is discussed in Appendix B; a full account of the waveform’s systematics in these injection experiments will be given elsewhere [In Prep.].

The recovery results with TaylorF2 are summarized in Tab. 1. Results from the threshold-mass analysis with maximum mass constrained to be larger than the mass of PSR J0348+0432 are also shown in Fig. 2 [See Appendix C for a similar plot without the maximum mass constraint]. The left panel shows for each injection the cumulative probability distribution to find $M>M_{\text{thr}}$; the vertical line marks the collapse threshold. The right panel shows for each injection the inferred mass divided by the inferred threshold mass $\left[M/M_{\text{thr}}\right]_{\rm rec}$ versus the same injected quantity, $\left[M/M_{\text{thr}}\right]_{\rm inj}$. Erroneous recoveries would populate the top-left and bottom-right quadrants of the plot. The plot shows that the inference returns the correct estimate of the prompt collapse for the majority of the injections. SFHo 1.46+1.46 is a borderline case for which $P_{\text{PC}}\sim 40\%$. However, we observe that a few simulations that led to a prompt collapse (denoted by solid lines and circles in the plots), were not recovered as such. This misclassification had occured already at the level of injection parameters, since their position on the x-axis (right panel) lies below unity, which is due to the inherent uncertainty on the estimation of $k_{\text{thr}}$ from fitting NR data. The resulting error is comparable in size to the statistical error of our parameter estimation.

Results from the threshold tidal parameter analysis are shown in Fig. 3. The left panel shows for each injection the probability that $\tilde{\Lambda}$ is smaller than a given value. The latter should be compared to Eq. (10), shown as a black solid line. The right panel summarizes the inference results in a way analogous to the Fig. 2. The threshold tidal parameter analysis incorrectly predicts few cases for BNS with $\tilde{\Lambda}\sim\tilde{\Lambda}_{\text{thr}}$. SFHo 1.40+1.40 ($\tilde{\Lambda}=334$) and ALF2 1.50+1.50 ($\tilde{\Lambda}=382$) are predicted with 61% and 29% probability of producing a NS remnant while the simulations indicate prompt BH formation. The two SLy-SOR binaries ($\tilde{\Lambda}\sim 401$) are predicted with $\sim 22$% probability of prompt collapse while the merger produces a HMNS. DD2 1.59+1.59 ($\tilde{\Lambda}=332$) is predicted to promplty form a BH with 66% probability, with the simulation resulting in a HMNS of very short life of ${\sim}3$ ms.

We conclude that both methods infer correctly the merger outcome of the simulations, except for few cases corresponding to BNS close to the collapse threshold. Excluding these cases (in which the answer given is anyway inconclusive), the mass threshold analysis with the maximum mass constraint better captures the formation of a NS remnant, while the threshold $\tilde{\Lambda}$ analysis captures more robustly the prompt collapse cases. The mass threshold analysis without the maximum mass constraint give instead results comparable to the threshold $\tilde{\Lambda}$ analysis. This can be understood by inspecting the EOS posteriors in Appendix D. The EOS inference from the inspiral data constrains more strongly the EOS at densities comparable to the maximum density of the individual NS in the binary. These densities are those that determine the individual masses and thus the $\Lambda_{A}$ parameters. Introducing a maximum mass constraint effectively corresponds to introduce a lower bound on the mass distribution of the individual NS (and on $k_{\text{thr}}$), thus lowering the collapse probability.

V Application to GW170817

We apply our analysis methods to data from the first detected BNS event GW170817, by postprocessing the publicly available posteriors of the LIGO-Virgo collaboration released with Abbott et al. (2018b, 2019). In all of the analysis set-ups, the NS spins are aligned with the orbital angular momentum and the spin magnitudes are restricted to the “low-spin” prior range $\chi_{1,2}\in[-0.05,0.05]$.

V.1 Results

For the threshold-mass method (Sec. III.1), we process the posteriors of the spectral parameters $(\gamma_{0},\gamma_{1},\gamma_{2},\gamma_{3})$ published in Abbott et al. (2018b),

• •

without imposing an implicit constraint on $M_{\text{TOV}}^{\text{max}}$;

• •

with the additional hard cut of $M_{\text{TOV}}^{\text{max}}\geq 1.97{\rm M_{\odot}}$, corresponding to a $1\sigma$ conservative estimate of the PSR J0348+0432 mass measurement $2.01\pm 0.04{\rm M_{\odot}}$ Antoniadis et al. (2013);

• •

with the additional probabilistic weight quantifying the probability of $M_{\text{TOV}}^{\text{max}}$ being heavier than the PSR J0348+0432 mass.

For the threshold-$\tilde{\Lambda}$ method, we process published posteriors on tidal parameters from a number of different analyses. In particular, we consider methods that extend the BBH parameter space by the matter-related parameters

• •

$(\tilde{\Lambda},\delta\tilde{\Lambda})$, using four different waveform models (IMRPhenomPv2NRtidal, IMRPhenomDNRtidal, SEOBNRT, TaylorF2)

• •

$\Lambda_{s}=(\Lambda_{1}+\Lambda_{2})/2$, the symmetric tidal parameter, using IMRPhenomPv2NRtidal, along with the use of the EOS-insensitive relation for the antisymmetric tidal parameter $\Lambda_{a}(\Lambda_{s},q)$ (see Abbott et al. (2018b)), which can then be mapped to $\tilde{\Lambda}$;

• •

$\vec{\gamma}=(\gamma_{0},\gamma_{1},\gamma_{2},\gamma_{3})$ parametrizing, the EOS, from which $\tilde{\Lambda}$ can be derived, with and without a constraint on $M_{\text{TOV}}^{\text{max}}$ (see parametrized EOS method of Abbott et al. (2018b)).

First, for the threshold-mass method we show in Figure 4 the joint posterior of total mass $M$ and the threshold mass $M_{\text{thr}}$ (left) as well as the cumulative distribution function of their ratio $M/M_{\text{thr}}$ (right), obtained with the theshold mass analysis with and without the constraint $M_{\text{TOV}}^{\text{max}}\geq 1.97$ ${\rm M_{\odot}}$. The latter plot should be interpreted as the probability of prompt collapse as a function of $k_{\text{thr}}$, if we pretented to be totally agnostic on $k_{\text{thr}}$. Without the maximum mass constraint, the collapse probability ranges from $P_{\text{PC}}\sim 0.2$ to $P_{\text{PC}}\sim 0.85$ for the expected range of $k_{\text{thr}}$ (see orange line and white region in the plot). Including the constraint $M_{\text{TOV}}^{\text{max}}\geq 1.97$ ${\rm M_{\odot}}$ strongly disfavours a prompt collapse: $P_{\text{PC}}=0$ if $k_{\text{thr}}>1.4$, growing up to $P_{\text{PC}}\sim 0.5$ if $k_{\text{thr}}\sim 1.3$, if for very soft EOS and NS compactness $C_{\text{max}}\sim 0.33$.

However, $k_{\text{thr}}$ is not an independent unknown parameter; using the results of Sec. II.1 we estimate the value of $k_{\text{thr}}$ and $M_{\text{thr}}$ from the EOS parameters $\vec{\gamma}$. The resulting posterior of $M/M_{\text{thr}}$ is plotted as a cumulative distribution in Fig. 5. Here too, we find a significant difference between the analyses with and without the $M_{\text{TOV}}^{\text{max}}$ constraint, that estimate the prompt collapse probability at $0.09$ and $0.59$ respectively. The reason is that the $M_{\text{TOV}}^{\text{max}}$ constraint removes part of the EOS parameter space that is too soft to support a NS mass of $1.97$ ${\rm M_{\odot}}$ (and will most likely predict a prompt collapse). The effect on $P_{\rm PC}$ is significant, since the recovered binary parameters of GW170817 happen to lie close to the prompt-collapse threshold.

Fig. 7 shows the prompt collapse probability obtained with the threshold-$\tilde{\Lambda}$ method (cf. Fig. 3, left panel) We find a prompt collapse probability between $P_{\text{PC}}\sim 43\%$ and $74\%$, depending on the waveform approximant used for the analysis and on the inference method employed. The data from the EOS inference employed also in the threshold-mass analysis give the smalleest prompt collapse probability as a result of imposing the maximum mass constraint. If the constraint is relaxed the probability grows to $69\%$. All the analysis performing inference on $\tilde{\Lambda}$ give prompt collapse probability between $54\%$ and $60\%$; the waveform approximants estimating the largest $\tilde{\Lambda}$ clearly give the smaller probabilities. The largest prompt collapse probability is obtained using the EOS-insensitive relations in the $\Lambda_{s}$ inference and employing the $\tilde{\Lambda}$ threshold. This can be understood as the combined effect of using IMRPhenomPv2NRtidal as a template waveform model, which tends to introduce a systematic bias favouring prompt collapse (see discussion in Appendix B) and not having a constraint on $M_{\text{TOV}}^{\text{max}}$.

Resulting values for the probability of prompt collapse $P_{\rm PC}$ from the above analyses are listed in Table 2. We observe that the parameters of GW170817 are measured around the threshold region both for the threshold-mass method and for the threshold-$\tilde{\Lambda}$ method. Thus, overall there seems to be no definite answer as to whether the BNS merger was followed by a prompt collapse to a BH. However, if we focus on the analyses where the $M_{\text{TOV}}^{\text{max}}$ constraint can be imposed, to account for the mass measurement of PSR J0348+0432, we see that the prompt-collapse hypothesis is strongly disfavoured.

We also point out that the GW170817 inference of tidal effects using various point-mass waveform approximants combined with NRTidal gives posteriors with a bimodal distribution peaked around $\tilde{\Lambda}\sim 200$ and $\tilde{\Lambda}\sim 600$ and support up to $\tilde{\Lambda}\sim 800$; while using TaylorF2 and EOB approximants it gives a single broader peak at $\tilde{\Lambda}\sim 300$ Abbott et al. (2019, 2018a). Independent analysis confirm these findings De et al. (2018); Dai et al. (2018); Radice and Dai (2019).

VI Conclusion

We proposed two methods to infer prompt black hole formation from the analysis of the inspiral gravitational wave of a binary neutron merger. Both methods rely on numerical-relativity models of the prompt collapse threshold for quasicircular and nonspinning binary neutron star merger. The methods are validated with a set of 17 injection and recovery experiments, and verified against data from numerical relativity simulations. All the signals were correctly recovered with the exception of few cases close or at the collapse threshold. Improving such cases will require more precise numerical relativity models and simulations. We conclude that our analysis could be robustly applied to GW170817-like signals (single events) captured by advanced LIGO-Virgo at designed sensitivity. We also point out that waveform systematics may introduce important biases in the near-threshold region.

Application of these two methods on the GW170817 data gives no definitive answer to whether the BNS merger was followed by a prompt collapse into a BH, as the recovered masses and tidal parameters of the binary appear to lie in the vicinity of the threshold. However, if a constraint is applied on the maximum irrotational NS mass supported by the EOS, that is compatible with the mass measurements of PSR J0348+0432 and PSR J0740+6620, then we observe a strong preference against the prompt collapse hypothesis.

Acknowledgements.

The authors thank the LIGO-Virgo matter group for discussions, Liang Dai and Kenta Hotokezaka for sharing some data and information. MB, SB, FZ acknowledge support by the EU H2020 under ERC Starting Grant, no. BinGraSp-714626. DR acknowledges support from a Frank and Peggy Taplin Membership at the Institute for Advanced Study and the Max-Planck/Princeton Center (MPPC) for Plasma Physics (NSF PHY-1804048). Data analysis for this paper was performed on the supercomputer ARCCA; we are grateful for computational resources provided by Cardiff University, and funded by STFC grant ST/I006285/1. Data analysis was performed on the Virgo “Tullio” server at Torino supported by INFN. Numerical relativity simulations were performed on the supercomputer SuperMUC at the LRZ Munich (Gauss project pn56zo), on supercomputer Marconi at CINECA (ISCRA-B project number HP10BMHFQQ); on the supercomputers Bridges, Comet, and Stampede (NSF XSEDE allocation TG-PHY160025); on NSF/NCSA Blue Waters (NSF AWD-1811236; on ARA cluster at Jena FSU.

Appendix A Prompt collapse threshold from NR data

We collect NR data for the prompt collapse threshold $k_{\text{thr}}$ from Hotokezaka et al. (2011); Bauswein et al. (2013); Zappa et al. (2018); Köppel et al. (2019) into Tab. 3.

In the first three references, the collapse threshold is estimated by performing simulations with a given EOS and different masses and then linearly interpolating between the two simulations that braket the threshold. The uncertainty on the $k_{\text{thr}}$ is thus determined not only by the grid resolution of the simulations but also by how close the two simulation bracket the threshold. Typical relative errors are at the level of 5% although, surprisingly, detailed grid resolution studies for this application are missing. In Köppel et al. (2019) the collapse threshold for a specific EOS is determined as the binary mass for which the merger remnant collapses over the free-fall timescale of a $M_{\text{TOV}}^{\text{max}}$ NS. The actual value of $M_{\text{thr}}$ is computed using an extrapolation of an exponential fit of the results obtained by a few prompt collapse simulations. The reported error is also the one obtained by the exponential fit. Despite providing consistent results, the error estimate is qualitatively different from the other approaches and counterintuitively the fewer the simulations the smaller the uncertainties. We thus select the $M_{\text{thr}}$ for which at least three simulations were performed and we assign to all of them the relative error obtained by their DD2 model, since this is the model with more points. We stress that this relative error is comparable to the smallest relative errors used in the other works.

The data are plotted as a function of the maximum compactness and shown in the top panel of Fig. 8. All the data show an approximate linear correlation with $C_{\text{max}}$ although part of the data of Hotokezaka et al. (2011) appear to systematically underestimate the threshold with respect to the other datasets. We perform several linear fits combining the datasets and propagating the uncertainties. The fit coefficients and their errors are reported in the bottom panel of Fig. 8. The different sets of coefficients are essentially compatible with each other, and the errors become smaller when more points are included. The equation we use in the main body (Eq. (3)) is the best fit given by the combination of the data of Bauswein et al. (2013); Zappa et al. (2018); Köppel et al. (2019).

Appendix B Waveform systematics

Our injection experiments highlight systematics biases in the recovery of the TEOBResumS waveforms when using the IMRPhenomPv2NRtidal as our template model, while results are overall consistent when TaylorF2 is employed with or without a high cut-off frequency of $1024$ Hz. Representive measurements of the tidal parameter $\tilde{\Lambda}$ are shown in Fig. 9 for two injections with different EOS. A similar bias between TEOBResumS and IMRPhenomPv2NRtidal is visible in the injections of Dudi et al. (2018); Messina et al. (2019) performed at SNR 100, but so far this has not been systematically investigated nor explained. We plan to do so in a forthcoming publication. Here, we use the IMRPhenomPv2NRtidal results to discuss how waveform systematics can affect the prompt collapse inference with our two methods. In the comparison plots of Fig. 9 we show how analyzing the same signal with different waveform models can affect the estimated prompt collapse probability; in particular IMRPhenomPv2NRtidal tends to underestimate tidal deformabilities and therefore overestimate $P_{\rm PC}$. A similar effect is observed in the threshold-mass results for the same injections in the left panel of Fig. 10, where $M/M_{\text{thr}}$ is overestimated by IMRPhenomPv2NRtidal (and thus so is $P_{\rm PC}$). Overall, we find the significance of waveform systematics to be limited to the cases close to the collapse threshold.

Figure 11 summarizes waveform systematics effects on the threshold tidal parameter analysis. We find that the prompt collapse inference with TaylorF2 and cut-off $1024$ Hz gives consistent results with the injection except for binaries with $\Lambda\sim\tilde{\Lambda}_{\text{thr}}$. For the DD2 1.59+1.59 BNS ($\tilde{\Lambda}=332$) and the two SLy binaries ($\tilde{\Lambda}\sim 401$) the method estimates respectively a 75% and $\sim 40$% probability of prompt collapse while the merger result in a HMNS. In the former case the binary is at the collapse threshold and the HMNS is very short lived ($3$ ms). Hence it could be simply a results of out uncertainties. In the latter case the binaries are slightly above the collapse threshold and the prediction appears to have a genuine systematic error of the method. Similarly, for SFHo 1.40+1.40 ($\tilde{\Lambda}=334$) and ALF2 1.50+1.50 ($\tilde{\Lambda}=382$) the method predicts 33% and 54% probability of producing a NS remnant while the simulations indicate prompt BH formation.

Recovering with IMRPhenomPv2NRtidal systematically underestimates the injected TEOBResumS $\tilde{\Lambda}$; the effect being worst for cut-off frequency $1024$ Hz and minimzed by cut-off $2048$ Hz. The result can be in part understood from the fact that the low frequency limit of the NRtidal is accurate only to the leading-order post-Newtonian tidal term Dietrich et al. (2017, 2019). The same systematic trend can be seen in the threshold-mass analysis summarized in the right panel of Fig. 10, which is more pronounced in the less compact binaries. The errors in the prompt collapse analysis due to the numerical fits on $k_{\text{thr}}$ discussed above, are now combined with those from the waveform systematics. As a result, the method predicts correctly the prompt collapse of ALF2 1.59+1.59 and SFHo 1.40+1.40 (thanks to a “cancellation” of systematic errors) but incorrectly favours prompt collapse for the SLy binaries.

Appendix C Effect of $M_{\text{TOV}}^{\text{max}}$ constraint

In the threshold-mass method, sampling the EOS parameter space directly allowed us to impose conditions on the maximum stable nonrotating NS mass, $M_{\text{TOV}}^{\text{max}}$. In this section we examine the effect that different choices of this constraint may have on estimating the probability of prompt collapse.

First we review the results of the injection study of Sec. IV when imposing a $M_{\text{TOV}}^{\text{max}}$ constraint based on the mass measurement of PSR J0348+0432. Results are summarized in Fig. 12. When comparing against Fig. 2, we observe a systematic trend to lower values of recovered $M/M_{\text{thr}}$. This can be interpreted as a push towards higher values of $M_{\text{thr}}$, which is expected, since a soft part of the space of EOS is effectively removed from our prior. Note the peculiar behavior of the 2B BNS as a consequence of the fact that the maximum mass for that EOS violates the prior imposed in the analysis.

We now move on to the analysis of GW170817 data using the spectral EOS parametrization and consider the following choices:

• •

No constraint on $M_{\text{TOV}}^{\text{max}}$;

• •

A hard constraint of $M_{\text{TOV}}^{\text{max}}\geq 1.97{\rm M_{\odot}}$, corresponding to a conservative 1-$\sigma$ bound on the mass of PSR J0348+0432;

• •

A probabilistic constraint based on the mass measurement of PSR J0348+0432, which follows the Gaussian PDF $\mathcal{N}(2.01,0.04)$;

• •

A probabilistic constraint based on the recent observation of PSR J0740+6620, which follows the Gaussian PDF $\mathcal{N}(2.17,0.11)$;

The results are illustrated in Fig. 13. We find that if the heavy-NS measurements are taken into account, the prompt-collapse probability tends to zero (even more so than in the case of a hard cut at $1.97\>{\rm M_{\odot}}$).

Appendix D EOS reconstruction

In our threshold-mass method we have employed the spectral family of Lindblom (2010) to parametrize the EOS. It is instructive to examine whether the posterior PDF on the EOS parameters $\vec{\gamma}$ faithfully reconstructs the injected model, within the margins of our measurement error. Two typical cases are illustrated in Fig. 14, where the reconstructed $P(\rho)$ curves are clearly distinguishable from each other and faithfully follow the corresponding underlying model. In particular we observe the separation becoming more clear around a “focal point” at $\sim 2\>\rho_{\rm nuc}$, which happens to be close to the typical central density of the NS, which largely determines the bulk properties of the star.