Using Gravitational-wave Observations and Quasi-universal Relations to Constrain the Maximum Mass of Neutron Stars

A long-awaited event took place on 2017 August 17: the Advanced LIGO and Virgo network of GW detectors have recorded the signal from the inspiral and merger of a binary neutron-star (BNS) system (Abbott et al. 2017b). The correlated electromagnetic signals that have been recorded by ∼70 astronomical observatories and satellites have provided the striking confirmation that such mergers can be associated directly with the observation of short gamma-ray bursts (SGRBs). This event has a double significance. First, it effectively marks the birth of multi-messenger GW astronomy. Second, it provides important clues to solve the long-standing puzzle of the origin of SGRBs (Eichler et al. 1989; Narayan et al. 1992; Rezzolla et al. 2011; Berger 2014). Numerical simulations in full general relativity of merging BNSs have also played an important role in determining the solution of this puzzle, and significant progress has been made over the last decade to accurately simulate the late-inspiral, merger, and post-merger dynamics of BNSs (see, e.g., Baiotti & Rezzolla 2017; Paschalidis 2017 for recent reviews).

Indeed, it is through the detailed analysis of the results of these simulations that a number of recent suggestions have been made on how to use the GW signal from merging BNSs to deduce the properties of the system and, in particular, the equation of state (EOS) of nuclear matter.

For instance, the changes in the phase evolution of the GW signal during the inspiral, which depends on the tidal deformability of stellar matter will leave a characteristic imprint on the GW signal (Read et al. 2013; Bernuzzi et al. 2014; Hinderer et al. 2016; Hotokezaka et al. 2016) or in the post-merger phase. This imprint, such as the one associated with the GW frequency at maximum amplitude (Read et al. 2013; Bernuzzi et al. 2014; Takami et al. 2015), can even be quasi-universal in the sense that it depends only weakly on the EOS. Similar considerations also apply for the post-merger signal, where the GW spectrum exhibits characteristic frequencies (Bauswein & Janka 2012; Takami et al. 2014), some of which have been shown to have a quasi-universal behavior (Bernuzzi et al. 2014; Takami et al. 2014, 2015; Rezzolla & Takami 2016; Maione et al. 2017).

Much more subtle, however, has been the task of determining the precise fate of the binary merger product (BMP), as this depends on a number of macroscopical factors, such as the total mass and mass ratio of the BNS system of the angular-velocity profile (Hanauske et al. 2017), but also of microphysical ones, such as the efficiency of energy transport via neutrinos (Palenzuela et al. 2015; Sekiguchi et al. 2016; Bovard et al. 2017) and the redistribution of angular momentum via magnetic fields (Siegel et al. 2014; Palenzuela et al. 2015; Endrizzi et al. 2016). While attempts have been made to determine the mass of the binary that would lead to a prompt collapse, i.e., to a black hole within few milliseconds after merger (see, e.g., Baiotti et al. 2008; Bauswein et al. 2013), or to determine the lifetime of the merged object (see, e.g., Lasky et al. 2014; Ravi & Lasky 2014; Piro et al. 2017), the picture on the fate of the post-merger object is still rather uncertain. What makes such a picture complicated is the multiplicity of stable, unstable, and metastable equilibria in which the merged object can find itself. The importance of clarifying this picture, however, is that understanding the ability of the merged object to sustain itself against gravitational collapse is directly related to the maximum mass that can be sustained against gravity, which depends on the underlying EOS.

In this Letter, we combine the recent GW observation of the merging system of BNSs via the event GW170817 (Abbott et al. 2017b) with the existence of quasi-universal relations regulating the equilibria of rotating and nonrotating compact stars to set constraints on the maximum mass that can be sustained by nonrotating stellar models of neutron stars. More specifically, after defining the maximum mass of nonrotating models, ${M}_{{}_{\mathrm{TOV}}}$, and recalling that the maximum mass that can be supported through uniform rotation is ${M}_{\max }=({1.20}_{-0.02}^{+0.02}){M}_{\mathrm{TOV}}$ independently of the EOS (Breu & Rezzolla 2016), we deduce that when the merged object collapses it has a core that is uniformly rotating and close to the maximum mass of uniformly rotating configurations. Then our range reduces considerably and sets the following constraint for the maximum mass ${2.01}_{-0.04}^{+0.04}\leqslant {M}_{\mathrm{TOV}}/{M}_{\odot }\lesssim {2.16}_{-0.15}^{+0.17}$. Our estimate, which is compatible with that recently suggested by other authors (Alsing et al. 2017; Margalit & Metzger 2017; Shibata et al. 2017; Ruiz et al. 2017), follows a straightforward set of considerations and does not rely on the modeling of the electromagnetic signal via numerical-relativity simulations (as done, e.g., by Bovard et al. 2017 or Shibata et al. 2017) but only on basic arguments inferred from kilonova modeling (Cowperthwaite et al. 2017), can be further refined as new observations are carried out.

To illustrate the multiplicity of states that the merger of a BNS system can lead to, we show in Figure 1 a schematic diagram reporting the (gravitational) mass ${M}_{{\rm{g}}}$ versus the central rest-mass density ${\rho }_{c}$ and thus comprising all possible stable and unstable equilibrium states for the BMP. More specifically, shown with two solid black lines are the sequences of nonrotating (bottom) neutron stars and the neutron stars spinning at the mass-shedding limit. The vertical thin black line marks the turning points of sequences with constant angular momentum and has been shown to be a good approximation to the neutral-stability line for uniformly as well as differentially rotating neutron stars (Takami et al. 2011; Weih et al. 2018). Models on the low-density side of this line are dynamically stable, while the ones on the high-density side are unstable against gravitational collapse to a black hole. Neutron stars with masses exceeding the maximum mass of nonrotating configurations, ${M}_{{}_{\mathrm{TOV}}}$, but not the maximum mass of uniformly rotating neutron stars, ${M}_{\max }$, are referred to as supramassive (SMNS), while the ones with mass higher than ${M}_{\max }$ are called hypermassive (HMNS; dark-red shaded area in Figure 1). The latter configurations can only be supported by differential rotation. SMNSs, on the other hand, can be either uniformly or differentially rotating. The uniformly rotating models, however, are confined to the region between the nonrotating and mass-shedding limit (green area). Outside this region, only differentially rotating SMNSs are possible (medium-red area). Finally, models with mass below ${M}_{\mathrm{TOV}}$ can be rotating either differentially (light-red area) or uniformly (light-green area).

Zoom In Zoom Out Reset image size

Also reported as dashed lines are two "trajectories" that the BMP produced in GW170817 (golden cross), could have followed and that we have indicated as $(1)$ and $(2)$, respectively. Both trajectories end on the neutral-stability line because we hereafter make the working assumption that the BMP produced in GW170817 has indeed collapsed to a black hole as is necessary for most models of SGRB emission from BNS mergers, see, e.g., Rezzolla et al. (2011) and Murguia-Berthier et al. (2017), and is also expected to occur for most commonly used EOSs given the total mass of the system (Abbott et al. 2017a).

In the first scenario $(1)$, the BMP spins down and redistributes its angular momentum due, for instance, to magnetic braking or the development of a magnetorotational instability (see Baiotti & Rezzolla 2017 for a review). It does so moving on a line of almost constant baryon mass until it eventually enters the dark-green shaded region on the stable side of the neutral-stability line. It can then further lose gravitational mass by spinning down until it eventually crossed the neutral-stability line as a uniformly rotating SMNS and collapses. In the second scenario $(2)$, instead, the BMP passes the neutral-stability line much more rapidly and before it can redistribute its angular momentum, thus collapses as a differentially rotating HMNS. This scenario, however, is unlikely when considering the blue-kilonova signal that has been observed (Cowperthwaite et al. 2017) in the electromagnetic counterpart of GW170817. To produce such a signal, in fact, ejected material with very high electron fraction ${Y}_{e}\gt 0.25$ must be produced, which, however, most likely originate from the hot polar region of the BMP (Bovard et al. 2017; Metzger 2017a, 2017b). Hence, the observation of such a signal inevitably requires the BMP to be sufficiently long-lived. In particular, its lifetime should be much longer than the timescale for reaching uniform rotation via magnetic braking.

These considerations make the scenario $(1)$ the most likely one. At the same time, the BMP cannot have survived for very long if an SGRB was observed only $\simeq 1\,{\rm{s}}$ after the merger, thus constraining the mass of the BMP to be very close to ${M}_{\max }$ when passing the neutral-stability line. This conclusion becomes inevitable when considering the timescales associated with the spinning down of a uniformly rotating neutron star. Magnetic-dipole emission, in fact, is not sufficiently efficient and would act on much longer timescales (see, e.g., Zhang & Mészáros 2001). Spin-down (and hence loss of gravitational mass) via the GW emission driven by an ellipticity in the BMP is of course possible, but would require unrealistic deformations to be efficient over only $1\,{\rm{s}}$. We reach this conclusion by estimating the ellipticity ε required to produce such a loss by considering the typical timescale of GW emission to be (Usov 1992)

$$\begin{eqnarray}&&{\tau }_{\mathrm{GW}}=\displaystyle \frac{{E}_{\mathrm{kin}}}{{L}_{\mathrm{GW}}},\end{eqnarray} \tag{ 1 }$$

where ${E}_{\mathrm{kin}}=I{{\rm{\Omega }}}^{2}/2$ and

$$\begin{eqnarray}&&{L}_{\mathrm{GW}}=\displaystyle \frac{32}{5}\displaystyle \frac{{{GI}}^{2}{{\rm{\Omega }}}^{6}}{{c}^{5}}{\varepsilon }^{2},\end{eqnarray} \tag{ 2 }$$

where I is the moment of inertia, Ω is the rotational frequency, G is Newton's constant, and c is the speed of light. Using typical values of $I\approx {10}^{45}\,{\rm{g}}\,{\mathrm{cm}}^{2}$, we find that

$$\begin{eqnarray}&&\varepsilon \gtrsim 3\times {10}^{-2}\ {\left(\displaystyle \frac{{10}^{4}{{\rm{s}}}^{-1}}{{\rm{\Omega }}}\right)}^{2}{\left(\displaystyle \frac{1{\rm{s}}}{{\tau }_{\mathrm{GW}}}\right)}^{\tfrac{1}{2}},\end{eqnarray} \tag{ 3 }$$

where we have intentionally underestimated the rotational frequency Ω of the remnant. Such high ellipticities are very unlikely even $30\,\mathrm{ms}$ after the merger since the BMP becomes essentially axisymmetric on timescales $\lesssim 50\,\mathrm{ms}$ (Hanauske et al. 2017). In summary, it is unlikely that the BMP has crossed the stability line as a differentially rotating object, as this would have happened on a timescale of tens of milliseconds. At the same time, it must have crossed the stable region for uniform rotation very rapidly, or it would have survived for timescales of the order of thousands of seconds. Hence, we conclude that it must have collapsed very close if not at the mass-shedding ${M}_{\max }$, which is what we will assume hereafter.

A way to exploit the information from GW observations to set constraints on the maximum mass of nonrotating stellar configurations (and hence on the EOS) has recently been suggested by the work of Breu & Rezzolla (2016). In that study, and inspired by the findings of Yagi & Yunes (2013), it was proposed that universal relations can be valid also away from regions of the space of stable solutions.

In particular, Breu & Rezzolla (2016) have shown that a universal relation is exhibited also by equilibrium solutions of rotating relativistic stars that are not stable. For this, uniformly rotating configurations on the turning-point line, i.e., whose mass is an extremum along a sequence of constant angular momentum, have been considered. Such configurations are unstable since they are found at larger central rest-mass densities than those on the neutral-stability line and are therefore marginally stable (Takami et al. 2011). In this way, it was possible to show that this relation holds not only for the maximum value of the angular momentum, but also for any rotation rate. The importance of this universal relation is that it allows one to compute the maximum mass sustainable through rapid uniform rotation, finding that, for any EOS, it is about 20% larger than the maximum mass supported by the corresponding nonrotating configuration, i.e., ${M}_{\max }\simeq ({1.20}_{-0.02}^{+0.02}){M}_{{}_{\mathrm{TOV}}}$, for all the EOSs considered. The existence of such a universal relation has been confirmed by several other authors and shown to apply also for other theories of gravity (see, e.g., Staykov et al. 2016; Minamitsuji & Silva 2016; Yagi & Yunes 2017).

Additionally, we show a quasi-universal relation for the conversion between gravitational mass and baryon mass, ${M}_{{\rm{b}}}$. In Figure 2, the conversion factor ${M}_{{\rm{b}}}/M$ is shown for the sequence of uniformly rotating neutron stars spinning at the mass-shedding limit. Interestingly, the value for the configuration with maximum mass is only weakly dependent on the underlying EOS, and we find

$$\begin{eqnarray}&&\eta := \displaystyle \frac{{M}_{{\rm{b}}}}{{M}_{\max }}\simeq 1.171,\end{eqnarray} \tag{ 4 }$$

with a standard deviation of $\sigma =6.8\times {10}^{-3}$. A similar universal relation for the conversion between baryon and gravitational mass has been proposed in Timmes et al. (1996) and Breu & Rezzolla (2016). We have used the estimate (4) here, as it refers specifically to models that are on the mass-shedding limit and have the maximum mass. This is presently the most accurate estimate possible for ${M}_{{\rm{b}}}/M$ at the mass-shedding limit and represents a considerable improvement over the relation derived by Timmes et al. (1996), i.e., ${M}_{{\rm{b}}}/M=1+0.075\,M$, which is shown as dashed lines in Figure 2. We note that this relation is often employed, e.g., by Piro et al. (2017), but it systematically overestimates the relation between the two masses by 10%–25%.

Zoom In Zoom Out Reset image size

We can now use these universal relations to derive a simple mass constraint on the EOS, making very basic assumptions on the mass distribution of the remnant.

As a simple parametrization, we assume that the system can be described by the amount of ejected baryon mass M_ej from the inner core of the remnant, the initial baryon mass M_b of the merger remnant and the baryon mass in the uniformly rotating core ${M}_{\mathrm{core}}=\xi {M}_{{\rm{b}}}$.

Now we can invoke simple baryon mass conservation to concluded that ${M}_{\mathrm{core}}(t=0)={M}_{\mathrm{core}}(t)+{M}_{\mathrm{ej}}$. As we have detailed in the previous sections, we assume that the remnant attains uniform rotation in the vicinity of the Keplerian limit, e.g., ${M}_{\mathrm{core}}\equiv {M}_{\mathrm{core}}({t}_{\mathrm{collapse}})={M}_{{\rm{b}},\max }$, where ${M}_{{\rm{b}},\max }=\eta {M}_{\max }$ is the baryon mass at the mass-shedding limit. Making the simplifying assumption and solving for M_max, we find

$$\begin{eqnarray}&&{M}_{\max }={\eta }^{-1}(\xi {M}_{{\rm{b}}}-\ {M}_{\mathrm{ej}}).\end{eqnarray} \tag{ 5 }$$

Combining this with the result from Breu & Rezzolla (2016), we infer

$$\begin{eqnarray}&&{M}_{\mathrm{TOV}}={\chi }^{-1}{M}_{\max }={\chi }^{-1}(\xi {M}_{{\rm{g}}}-{\eta }^{-1}\ {M}_{\mathrm{ej}}),\end{eqnarray} \tag{ 6 }$$

where $\chi ={1.20}_{-0.02}^{+0.02}$ (Breu & Rezzolla 2016) and ${M}_{{\rm{g}}}={\eta }^{-1}{M}_{{\rm{b}}}={2.74}_{-0.01}^{+0.04}$, which is consistent with low-spin priors (Abbott et al. 2017b).

The assumption that the core collapses exactly at the maximum mass-shedding limit, i.e., $\chi \simeq 1.2$, brings in an error that needs to be accounted for, by considering a lower value for χ (Equation (12) in Breu & Rezzolla 2016). We thus set the lower bound to $\chi =1.15$, corresponding to a star close to, but not at the maximum mass-shedding limit.

Hanauske et al. (2017) have found that the mass fraction of the core after dynamical mass ejection is roughly $\xi ={0.95}_{-0.06}^{+0.06}$ (see Table II in Hanauske et al. 2017). The mass of the ejecta from the core is harder to estimate but, using standard kilonova models (Metzger 2017b; Shibata et al. 2017), it is reasonable to associate them with the blue ejecta ${M}_{\mathrm{ej}}^{\mathrm{blue}}={0.014}_{-0.010}^{+0.010}$ (Cowperthwaite et al. 2017; Drout et al. 2017), where we have assumed a conservative kilonova model dependent error that we use as $2\sigma$ for assigning a Gaussian probability distribution to the blue ejecta.

The resulting fit for M_TOV is shown in Figure 3, where the dashed lines refer to errors in ξ and the red shaded region is modeled with a Gaussian distribution taking into account the errors of M_ej. This region is framed by the 90% credibility levels of the binary mass (Abbott et al. 2017b).

Zoom In Zoom Out Reset image size

In summary, collecting all available information, we conclude that the maximum mass that can be supported against gravity by a compact nonrotating star is in the range

$$\begin{eqnarray}&&{2.01}_{+0.04}^{+0.04}\lt {M}_{\mathrm{TOV}}/{M}_{\odot }\lt {2.16}_{-0.15}^{+0.17},\end{eqnarray} \tag{ 7 }$$

where the lower limit in the range (7) is actually derived from accurate observations of massive pulsars in binary systems (Antoniadis et al. 2013).

The error corresponds to twice the standard deviation ($\sim 90 \%$ confidence) computed with standard error propagation, where the asymmetric errors in M_g and χ are taken into account by computing the standard deviation for the upper and lower limit separately. Clearly, values close to the upper and lower limits are unlikely, given the fact that not all the values of M_g and M_ej are equally likely (compare to the red shaded area).

Note the interesting general trend shown by the maximum mass in Figure 3: the estimates for M_TOV grow systematically with increasingly massive binary systems and with decreasing ejected masses (compare to the shading from light to dark blue). Hence, future detections of merging binary systems with masses smaller than that of GW170817 will help set even tighter constraints on the maximum mass M_TOV.

We have combined the recent GW observations of merging systems of binary neutron stars via the event GW170817 with a quasi-universal relation between the maximum mass of nonrotating stellar models ${M}_{\mathrm{TOV}}$ and the maximum mass that can be supported through uniform rotation to set new and tighter constraints on ${M}_{\mathrm{TOV}}$.

Our estimate follows a simple line of arguments and is based on a single and reasonable assumption that the product of the merger measured with GW170817 has collapsed to a rotating black hole when it had reached a mass close to the maximum mass for SMNS models. In this way, we can exploit quasi-universal relations to deduce that the maximum mass for nonrotating stellar configurations should be in the range $2.012.{01}_{-0.04}^{+0.04}\leqslant {M}_{\mathrm{TOV}}/{M}_{\odot }\lesssim {2.16}_{-0.15}^{+0.17}$. We note that it is, in principle, possible to constrain the lower limit for ${M}_{\mathrm{TOV}}$ also with a quasi-universal relation on the maximum mass of a neutron star in differential rotation (Weih et al. 2018).

A few remarks before concluding. First, a much more conservative upper limit ${M}_{\mathrm{TOV}}$ can be set uniquely assuming that the maximum nonrotating mass ${M}_{\mathrm{TOV}}$ cannot be smaller than the mass in the uniformly rotating core M_core. Taking into account the amount of mass ejected and the conversion between baryon and gravitational mass, this yields ${M}_{\mathrm{TOV}}/{M}_{\odot }\lesssim 2.59$. Second, our predictions are compatible with those recently presented by Shibata et al. (2017) and Margalit & Metzger (2017), sharing a number of similar considerations with the latter. However, differently from these other works, we have not employed a simple correlation between the maximum mass-shedding mass and the maximum nonrotating mass, or fitting formulas stemming from numerical simulations whose error budget is uncertain (Bauswein et al. 2013), nor have we relied on direct comparisons with numerical-relativity simulations for the electromagnetic emission. Rather, using basic arguments from kilonova modeling (Cowperthwaite et al. 2017), we have exploited the power of universal relations for the maximum mass that are valid for any value of the specific angular momentum (Breu & Rezzolla 2016). Third, the results presented here already have a direct impact on some of the EOSs describing matter at nuclear densities (see, e.g., Oertel et al. 2017 for a recent review). For instance, a popular EOS routinely employed in numerical-relativity calculations such as the DD2 EOS (Typel et al. 2010), violates the constraint (7) since it has ${M}_{\mathrm{TOV}}=2.419\,{M}_{\odot };$ at the same time, EOSs with hyperons, e.g., BHBΛΦ (Banik et al. 2014) and DD2Y (Marques et al. 2017), have maximum masses $\lesssim 2.1\,{M}_{\odot }$ and therefore seem favored (Richers et al. 2017). Finally, we note that the procedure outlined here and the use of stacking techniques, as those developed in the analysis of the GW signal of BNSs (Del Pozzo et al. 2013; Agathos et al. 2015; Clark et al. 2016; Bose et al. 2017), can be employed in the future as the results of new detections become available to set new and tighter constraints on the maximum mass. New observations, in fact, will set sharper boundaries in the probability distributions presented in Figure 3, thus tightening the estimates for the maximum mass.

It is a pleasure to thank the referee for useful suggestions and Luke Bovard and Enping Zhou for discussions. Support comes in part from "NewCompStar," COST Action MP1304; LOEWE-Program in HIC for FAIR; European Union's Horizon 2020 Research and Innovation Programme (Grant 671698) (call FETHPC-1-2014, project ExaHyPE), the ERC Synergy Grant "BlackHoleCam: Imaging the Event Horizon of Black Holes" (Grant No. 610058).