A General-relativistic Determination of the Threshold Mass to Prompt Collapse in Binary Neutron Star Mergers

The recent detection of gravitational waves from the merger of neutron-star binaries (Abbott et al. 2017) has heralded the new era of multi-messenger gravitational-wave astronomy. These observations offer new insight into the most extreme objects in universe, namely, neutron stars, and allow us to probe and constrain the properties of nuclear matter (Bauswein et al. 2017; Paschalidis et al. 2018; Annala et al. 2018; Burgio et al. 2018; Most et al. 2018; Montana et al. 2018; Radice et al. 2018).

When two neutron stars merge, they will produce an object that either collapses promptly to a black hole, or does not (Baiotti et al. 2008). In the latter case, the remnant may be a metastable object, e.g., a hypermassive neutron star (HMNS), eventually collapsing to a black hole on a secular timescale, or survive for much longer times, either as a rotating or a nonrotating star (see, e.g., Baiotti & Rezzolla 2017, for a review). In the case of the first detection of merging neutron stars, GW170817 (Abbott et al. 2017), the precise fate of the merger remnant is presently unknown, although the formation of a black hole naturally matches the simultaneous observation of a short gamma-ray burst (Eichler et al. 1989; Rezzolla et al. 2011), and has been the working hypothesis to set new limits on the maximum mass of neutron stars (Margalit & Metzger 2017; Shibata et al. 2017; Rezzolla et al. 2018; Ruiz et al. 2018).

Determining the time of collapse of the merger remnant is particularly challenging as there are a number of physical processes that either determine or undermine the stability of merger remnant. These include: the ejection of matter (Rosswog et al. 1999; Kyutoku et al. 2014; Lehner et al. 2016; Radice et al. 2016; Bovard et al. 2017; Dietrich & Ujevic 2017), the angular-momentum transfer via magnetic fields (Siegel et al. 2013; Kiuchi et al. 2015; Kawamura et al. 2016), the evolution of the degree of differential rotation (Kastaun et al. 2016; Hanauske et al. 2017), and possible viscous effects mediated either by neutrinos or magnetic fields (Duez et al. 2004; Radice 2017; Shibata & Kiuchi 2017; Alford et al. 2018).

The determination of the critical (threshold) mass to a prompt collapse, M_th, is much simpler, although it still poses numerical and conceptual challenges. Bauswein et al. (2013) were the first to explore this problem by employing a smooth-particle approximation for the hydrodynamics and a conformally flat approximation to general relativity. In this way, they were able to find a linear universal relationship between M_th and the compactness of the maximum-mass model, ${{ \mathcal C }}_{{}_{\mathrm{TOV}}}:= {M}_{{}_{\mathrm{TOV}}}/{R}_{{}_{\mathrm{TOV}}}$, where ${M}_{{}_{\mathrm{TOV}}}$ and ${R}_{{}_{\mathrm{TOV}}}$ are, respectively, the mass and radius of the maximum-mass nonrotating star. Here, we improve on this result by using a fully general-relativistic approach, a wider range of compactnesses, and a rigorous definition of the threshold mass. As a result, we find a nonlinear relation between M_th and ${{ \mathcal C }}_{{}_{\mathrm{TOV}}}$, which offers a better match to the numerical-relativity results. Furthermore, exploiting the information from GW170817, we use the new relation to set more stringent lower bounds on the radii neutron stars.

To describe the evolution of the merging system, we solve the coupled Einstein-hydrodynamics system (Rezzolla & Zanotti 2013) using the Einstein Toolkit (Löffler et al. 2012). In particular, we evolve the spacetime with the McLachlan code (Brown et al. 2009), with the same gauges as in Hanauske et al. (2017). On the other hand, we evolve the matter with the high-order relativistic-hydrodynamics code WhiskyTHC (Radice & Rezzolla 2012; Radice et al. 2014). The numerical grid uses the fixed-mesh refinement driver Carpet (Schnetter et al. 2004), with a total of six refinement levels having a highest resolution of 215 m covering the two stars and a total extent of 700 km. For one equation of state (EOS), we also considered different resolutions of 215, 287 and 573 m, obtaining threshold masses within a variance of ${\rm{\Delta }}{M}_{\mathrm{th}}\lesssim 0.005\,{M}_{\odot }$.

A number of zero-temperature ("cold") EOSs are available for numerical simulations and rely on nuclear physics calculations following a variety of theoretical approaches. While these EOSs are suitable to describe the inspiral phase, they become obviously inadequate after the merger, when the temperatures reach values of several tens of MeV. To counter this, it is customary to model the post-merger dynamics by modifying these EOS and adding a "thermal" contribution via an ideal-fluid EOS (Rezzolla & Zanotti 2013) so as to account for the shock heating (Janka et al. 1993). This approach, while not self-consistent, is rather robust and such "hybrid EOSs" have been employed extensively in the literature (Baiotti & Rezzolla 2017).

However, because our goal here is to determine as accurately as possible the threshold mass to prompt gravitational collapse, it is essential that the description of the thermal effects in the matter is as realistic and self-consistent as possible. In turn, this forces us to consider EOSs that have a physically consistent dependence on temperature. Unfortunately, the number of EOSs that can be employed for this scope and that do not violate some basic nuclear physics requirement (as it is the case for the widely employed LS220 EOS; Tews et al. 2017), is much more restricted. Here, we have employed all of the five "hot" EOSs that can be used with confidence to determine the threshold mass. The corresponding properties, when expressed in terms of the masses and radii of the maximum-mass nonrotating configuration (hereafter TOV), are reported in Table 1. Similarly, Figure 1 shows the masses and radii of the TOV equilibrium solutions, both stable (solid lines) and unstable (dashed lines). Solid squares mark the maximum-mass solutions, while open circles represent the models used.

Zoom In Zoom Out Reset image size

Table 1.  Properties of the Maximum-mass Models for the EOSs Considered Here: BHBΛΦ (Banik et al. 2014), (HS-)DD2 (Typel et al. 2010), SFHo (Steiner et al. 2013), (HS-)TM1 (Hempel et al. 2012), and Togashi+ (Togashi et al. 2016)

EOS	${M}_{{}_{\mathrm{TOV}}}$	${R}_{{}_{\mathrm{TOV}}}$	${C}_{\mathrm{TOV}}$	${\tau }_{{}_{\mathrm{TOV}}}$	Mth	ΔMth
	(M⊙)	(km)		(μs)	$({M}_{{}_{\mathrm{TOV}}})$	$({M}_{{}_{\mathrm{TOV}}})$
BHBΛΦ	2.10	11.64	0.26	83.31	1.503	0.005
DD2	2.42	11.94	0.30	80.60	1.364	0.020
SFHo	2.06	10.34	0.29	70.44	1.391	0.016
TM1	2.22	12.6	0.26	91.70	1.520	0.015
Togashi+	2.23	10.17	0.32	66.12	1.298	0.000

Note. Reported are: the maximum mass of a nonrotating Star ${M}_{{}_{\mathrm{TOV}}}$, the corresponding radius ${R}_{{}_{\mathrm{TOV}}}$, the compactness ${C}_{{}_{\mathrm{TOV}}}$, and the freefall-timescale ${\tau }_{{}_{\mathrm{TOV}}}$. Also shown are the threshold masses M_th and the corresponding errors.

Download table as:  ASCIITypeset image

Using these EOSs, we have modeled the initial data under the assumption of irrotational quasi-circular equilibrium (Gourgoulhon et al. 2001) and computed it via the LORENE library for a total of 15 equal-mass (i.e., q = 1) binaries. The initial separation is of 45 km, so that the binaries perform around five orbits before the merger. We note that because the threshold mass for equal-mass binaries is always larger than for unequal-mass binaries, i.e., M_th(q = 1) > M_th(q < 1), the use of equal-mass binaries is not a restriction but rather optimizes the search for M_th (see also Bauswein et al. 2017, for a discussion).

The measurement of the threshold mass inevitably imposes a clear definition, but also a straightforward procedure to extract this information from the numerical-relativity simulations. Quite generically, one expects that the lifetime of the HMNS, or collapse time ${t}_{\mathrm{coll}}$, will decrease as the mass of the binary is increased, so that the threshold mass will represent the shortest possible lifetime. This definition, is however inconvenient as different EOSs will yield different "shortest lifetimes" and comparing different EOSs may introduce clear biases. We resolve this problem by building our analysis around two important logical steps. First, we consider the collapse time as a dimensionless quantity by expressing it in terms of the freefall time τ_ff that, for an object of mass M and radius R, is given by (Rezzolla & Zanotti 2013)

$$\begin{eqnarray}&&{\tau }_{\mathrm{ff}}(M,R):= \displaystyle \frac{\pi }{2}\sqrt{\displaystyle \frac{{R}^{3}}{2M}}.\end{eqnarray} \tag{ 1 }$$

Because R ∼ 1/M for stable models, the smallest freefall time will be achieved for the maximum-mass model, so that the shortest freefall timescale is ${\tau }_{{}_{\mathrm{TOV}}}:= {\tau }_{\mathrm{ff}}({M}_{{}_{\mathrm{TOV}}},{R}_{{}_{\mathrm{TOV}}})$. Second, we define the threshold mass M_th as the one for which the merger remnant will collapse over such a timescale, i.e., $M/{M}_{{}_{\mathrm{TOV}}}\to {M}_{\mathrm{th}}/{M}_{{}_{\mathrm{TOV}}}$ for ${t}_{\mathrm{coll}}/{\tau }_{{}_{\mathrm{TOV}}}\to 1$.

We next discuss the procedure followed for the measurement of the collapse time, t_coll as given by the difference between the (coordinate) time of merger, t_merg, and that of collapse to a black hole, ${t}_{{}_{\mathrm{BH}}}$, i.e., ${t}_{\mathrm{coll}}:= {t}_{{}_{\mathrm{BH}}}-{t}_{\mathrm{merg}}$. Both of these times can be measured in a number of different ways. The first one involves the emission of gravitational waves, with t_merg being given by the time of the first maximum of the gravitational-wave strain amplitude h, or of the Weyl scalar Ψ₄. The time of collapse, on the other hand, can be estimated as the time when the ringdown signal starts. Overall robust measurements involving gravitational waves are prone to errors as t_merg depends sensitively on the phase evolution of the binary in its most nonlinear stage.¹ Similarly, determining ${t}_{{}_{\mathrm{BH}}}$ is complicated by the fact that the beginning of the ringdown is somewhat arbitrary and ringdown itself can be modified by the matter infalling onto the black hole.

To counter these difficulties, t_merg could be measured via the proper separation between the two stellar cores, marking the merger as the time when such a separation is below a fraction of the initial diameters of the two stars. Similarly, ${t}_{{}_{\mathrm{BH}}}$ could also be measured in terms of the first appearance of an apparent horizon. Although these measurements show an overall consistency, they suffer from the fact that the proper separation is very sensitive to the properties of the EOSs, and that the first appearance of the apparent horizon is ultimately set by the frequency at which it is searched during the simulations.

A third approach for measuring t_coll and ${t}_{{}_{\mathrm{BH}}}$ instead involves the monitoring of the minimum of the lapse function α, which we evolve employing a singularity-avoiding "$1+\mathrm{log}$" slicing (Alcubierre 2008). This quantity has been shown to be a very good proxy for the tracking and appearance of an apparent horizon (Alcubierre 2008) and has the advantage of being extremely robust. As a result, we can mark the two times respectively when

$$\begin{eqnarray}\begin{array}{rcl}{t}_{\mathrm{merg}}:\min (\alpha ) & = & {\alpha }_{\mathrm{merg}}:= 0.35,\\ {t}_{{}_{\mathrm{BH}}}:\min (\alpha ) & = & {\alpha }_{{}_{\mathrm{BH}}}:= 0.2.\end{array}\end{eqnarray} \tag{ 2 }$$

These values for $\min (\alpha )$ are less arbitrary than they may appear as α_merg =  0.35 and ${\alpha }_{{}_{\mathrm{BH}}}=0.2$ systematically represent the first minimum and zero of second derivative of the function $\min (\alpha (t))$, respectively; also, at ${\alpha }_{{}_{\mathrm{BH}}}=0.2$ the apparent horizon is also normally first found. More importantly, the results are invariant under a change of the values in (2).

In summary, all methods to compute t_coll provide results that are consistent across different choices. However, when it comes to robustness and simplicity of implementation, the monitoring of the minimum of the lapse function represents the optimal choice and is the one employed in the results that will be presented next.

Figure 2 reports the collapse times normalized to the freefall timescale of the maximum-mass models, ${t}_{\mathrm{coll}}/{\tau }_{{}_{\mathrm{TOV}}}$, for all of the EOSs considered; these times are then shown as a function of the initial half-mass of the binary normalized to the maximum-mass, M/M_TOV. The adoption of such set of dimensionless quantities has the goal of revealing a behavior of the threshold mass that is universal, i.e., only weakly dependent on the EOS (see, e.g., Breu & Rezzolla 2016; Weih et al. 2018, for some examples).

Zoom In Zoom Out Reset image size

Filled circles of different colors in Figure 2 report the numerical data for the various EOSs and mark the values of the initial masses in the binaries. Note that as t_coll decreases, even small differences in the initial masses can lead to rather large differences in the survival time. Also note that only two values are reported for the Togashi+ EOS, and these differ by only 3.7% in mass (i.e., M = 1.440, 1.435 M_⊙); a binary with a slightly smaller mass (e.g., M = 1.430 M_⊙) leads to an HMNS effectively stable over the timescales investigated here. Finally, as t_coll should diverge for vanishingly small values of M, we fit the numerical data with a simple exponentially decaying function of the type $M/{M}_{{}_{\mathrm{TOV}}}=\tilde{a}\exp [-\tilde{b}{({t}_{\mathrm{coll}}/{\tau }_{{}_{\mathrm{TOV}}})}^{2}]$. A rapid inspection of Figure 2 shows that the exponential fit is rather good and that a linear approximation would overestimate the threshold mass. Indeed, near the freefall limit, the behavior of $M/{M}_{{}_{\mathrm{TOV}}}$ should not be linear as in this limit infinitesimal changes in M are sufficient enough to yield a prompt collapse. Such a behavior, frequently encountered in critical-collapse calculations, requires the function $M/{M}_{{}_{\mathrm{TOV}}}$ to have vanishing derivative for ${t}_{\mathrm{coll}}/{\tau }_{{}_{\mathrm{TOV}}}\to 1$. Clearly, our nonlinear fitting reflects this behavior, while a linear one does not.

Figure 2 also reveals that the threshold mass is roughly given by

$$\begin{eqnarray}&&\displaystyle \frac{{M}_{\mathrm{th}}}{{M}_{{}_{\mathrm{TOV}}}}\approx 1.415,\end{eqnarray} \tag{ 3 }$$

with an uncertainty of ΔM_th = 0.05 M_⊙, i.e., a relative error of ∼4%. Hence, Equation (3) provides the lowest-order approximation between the threshold mass and the corresponding maximum mass. The existence of a relation of this type has been suggested initially by Bauswein et al. (2013), who concluded that ${M}_{\mathrm{th}}/{M}_{{}_{\mathrm{TOV}}}=k$, with k a linear function of the maximum compactness ${{ \mathcal C }}_{{}_{\mathrm{TOV}}}$. More specifically, the universal ansatz proposed by Bauswein et al. (2013) is that $k=\hat{a}\,{{ \mathcal C }}_{{}_{\mathrm{TOV}}}+\hat{b}$, with $\hat{a}=3.38,\hat{b}=2.43$, independent of the EOS.

Such a linear ansatz does represent a reasonable first approximation to the data, but not the most general one. In particular, if it represents a universal behavior of compact self-gravitating objects it should be valid for all possible compactnesses and provide the expected black hole limit, for which ${M}_{\mathrm{th}}/{M}_{{}_{\mathrm{TOV}}}\to 0$ for ${{ \mathcal C }}_{{}_{\mathrm{TOV}}}\to 1/2$. Hence, we correct the linear approximation via a nonlinear fit of the type

$$\begin{eqnarray}&&\displaystyle \frac{{M}_{\mathrm{th}}}{{M}_{{}_{\mathrm{TOV}}}}=a-\displaystyle \frac{b}{1-c\,{{ \mathcal C }}_{{}_{\mathrm{TOV}}}},\end{eqnarray} \tag{ 4 }$$

where a, b, c are to be determined from the data. However, imposing the fulfillment of the black hole limit removes one free parameter and sets $a=2b/(2-c)$.

Figure 3 reports with a solid blue line the fit (4) with b = 1.01, c = 1.34, against the numerical-relativity data shown with stars of the same colorcode as in Figure 2 (see Table 1 for the errors in the fit). Also shown with a red dashed line is the linear approximation of Bauswein et al. (2017), which clearly suggests larger threshold masses, most probably because the conformally flat approximation used by Bauswein et al. (2013) underestimates the strong-curvature behavior that characterizes the threshold to black hole collapse. At the same time, the difference with the linear approximation of is not enormous, and is 8% at most for the cases considered here.

Zoom In Zoom Out Reset image size

Additionally, we show in the inset of Figure 3 a comparison of the linear and nonlinear fittings for the threshold mass in the whole range of possible compactnesses, i.e., ${{ \mathcal C }}_{{}_{\mathrm{TOV}}}\in [0,1/2]$. The light green area refers to neutron stars, with the lower limit ${{ \mathcal C }}_{{}_{\mathrm{TOV}}}\gtrsim 0.2$ deduced from a large sample of EOSs (Most et al. 2018) and the upper limit ${{ \mathcal C }}_{{}_{\mathrm{TOV}}}\lesssim 0.35$ set by causality; instead, the upper limit of the dark green area ${{ \mathcal C }}_{{}_{\mathrm{TOV}}}\lt 4/9$ is set by the Buchdahl limit for compact stars (Rezzolla & Zanotti 2013).

In addition to providing an improvement over the linear approximation of Bauswein et al. (2013), the nonlinear expression (4) can now be used to provide more stringent (i.e., larger) lower limits on the radii of possible stellar models in the light of the recent detection of GW170817 (Abbott et al. 2017). In particular, following Bauswein et al. (2017), we use (4) to plot the threshold mass for different potential choices of ${R}_{{}_{\mathrm{TOV}}};$; this is shown in the left panel of Figure 4 with black solid lines and for ${R}_{{}_{\mathrm{TOV}}}=10,11,12\,\mathrm{km}$. Also reported in Figure 4 with a gray-shaded area is the limit set by causality and that requires ${M}_{{}_{\mathrm{TOV}}}/{R}_{{}_{\mathrm{TOV}}}\lesssim 0.354$ (Koranda et al. 1997). As noted by Bauswein et al. (2017), given the merger of a neutron-star binary with total mass M_tot, it is possible to set a lower limit on ${M}_{\mathrm{th}}$. This is shown in the left panel Figure 4, where we report with a horizontal blue dashed line the total gravitational mass estimated for GW170817, ${M}_{\mathrm{tot}}={2.74}_{-0.01}^{+0.04}\,{M}_{\odot }$ (Abbott et al. 2017).

Zoom In Zoom Out Reset image size

The corresponding uncertainty band (blue shaded area) gives a lower constraint on M_th, because GW170817 did not lead to a prompt collapse. The blue band thus constraints the red-shaded area from below, yielding a lower limit for the radius of the maximum-mass star, ${R}_{{}_{\mathrm{TOV}}}\geqslant {9.74}_{-0.04}^{+0.14}\,\mathrm{km}$ (red solid line); this is to be contrasted with the value deduced by Bauswein et al. (2017), i.e., ${R}_{{}_{\mathrm{TOV}}}\geqslant {9.26}_{-0.03}^{+0.17}\,\mathrm{km}$, on the basis of their linear approximation. Interestingly, to obtain a more stringent constraint similar to the one derived here, Bauswein et al. (2017) required a hypothetical detection of a binary with ${M}_{\mathrm{tot}}\simeq 2.9\,{M}_{\odot }$.

The logical approach and the mathematical procedure followed so far to derive the nonlinear fit (4) for ${{ \mathcal C }}_{{}_{\mathrm{TOV}}}$ can be repeated for the compactness of a fixed mass M_x, i.e., ${{ \mathcal C }}_{x}:= {M}_{x}/{R}_{x}$, thus allowing us to set constraints not only on ${R}_{{}_{\mathrm{TOV}}}$, but on any radius R_x within a reasonable range (the fit becomes increasingly bad for large masses as the EOSs do not have solutions well above two solar masses). This is shown in the right panel of Figure 4, where the values of M_x and R_x are indicated with blue crosses, while in black is the quadratic fit

$$\begin{eqnarray}&&{R}_{x}=-0.88\,{M}^{2}+2.66\,M+8.91.\end{eqnarray} \tag{ 5 }$$

The importance of expression (5) is that it offers a handy expression for the lower limit of stellar models as deduced from GW170817. A similar procedure has been followed by Bauswein et al. (2017), but only for a fixed mass of 1.6 M_⊙, from which it was deduced that R_1.6 ≥ 10.30 km; this result should be contrasted with the value derived from (5), which is instead R_1.6 ≥ 10.90 km. Similarly, for a reference star of 1.4 M_⊙ we obtain R_1.4 ≥ 10.92 km, which is worryingly close to the estimate by Bauswein et al. (2017) for 1.6 M_⊙.² On the other hand, our R_1.4 estimate is in good agreement with that of Most et al. (2018), who have explored a large number of possible EOSs and built a set of one billion stellar models from which they deduced that 12.00 < R_1.4/km < 13.45.

The detection of the merger of a binary system of neutron stars has concretely initiated the process of extracting information on the most extreme state of matter from gravitational-wave signals. An important phase of this process lies in understanding the post-merger behavior of the binary and the stability of the remnant. Indeed, knowing whether the system promptly forms a black hole is critical to understanding and interpreting the electromagnetic signals that may be observed.

Using a fully general-relativistic approach and a novel method for the determination of the threshold mass, we have carried out simulations making use of all of the realistic EOSs available to describe this process. In this way, we have found a nonlinear universal relation for the threshold mass as a function of the maximum compactness that is potentially valid for all compactnesses. At least for the temperature-dependent EOSs considered here, this universal relation improves the linear relation found recently with methods that are less accurate, but that also yield quantitatively similar results. Furthermore, exploiting the detection of GW170817, we have used the universal relation to set lower limits on the stellar radii for any any mass.

These results can be improved in at least two ways. First, as new hot EOSs become available for numerical simulations it will be possible to extend the analysis carried here, reducing its uncertainty. Second, as new detections from binary neutron-star mergers are revealed, the masses of these systems and their electromagnetic counterparts will be used to set ever more precise lower bounds on the radii of neutron stars.

Support comes from the ERC synergy grant "BlackHoleCam: Imaging the Event Horizon of Black Holes" (grant No. 610058), "PHAROS," COST Action CA16214; the LOEWE-Program in HIC for FAIR; the European Union's Horizon 2020 Research and Innovation Programme (grant No. 671698) (call FETHPC-1-2014, project ExaHyPE). The simulations were performed on the clusters SuperMUC (LRZ, Garching), LOEWE (CSC, Frankfurt), and HazelHen (HLRS, Stuttgart). The EOSs employed can be found on stellarcollapse.org.