Constraining the Maximum Mass of Neutron Stars from Multi-messenger Observations of GW170817

On 2017 August 17, the Advanced Laser Interferometer Gravitational-wave Observatory (LIGO) and Virgo network of gravitational-wave (GW) observatories discovered the inspiral and coalescence of a binary neutron star (BNS) system (LIGO Scientific Collaboration & Virgo Collaboration 2017), dubbed GW170817. The measured binary chirp mass was ${{ \mathcal M }}_{{\rm{c}}}={1.118}_{-0.002}^{+0.004}\,{M}_{\odot }$, with larger uncertainties on the mass of the individual neutron star (NS) components and total mass of ${M}_{1}=1.36\mbox{--}1.60\,{M}_{\odot }$, ${M}_{2}=1.17\mbox{--}1.36\,{M}_{\odot }$, and ${M}_{\mathrm{tot}}={M}_{1}\,+{M}_{2}={2.74}_{-0.01}^{+0.04}\,{M}_{\odot }$, respectively. These masses are derived under the prior of low dimensionless NS spin ($\chi \lesssim 0.05$), characteristic of Galactic BNS systems.

The electromagnetic follow-up of GW170817 was summarized in LIGO Scientific Collaboration et al. (2017a). The Fermi and INTEGRAL satellites discovered a sub-luminous gamma-ray burst (GRB) with a sky position and temporal coincidence within ≲2 s of the inferred coalescence time of GW170817 (Goldstein et al. 2017; LIGO Scientific Collaboration et al. 2017b; Savchenko et al. 2017). Eleven hours later, an optical counterpart was discovered (Allam et al. 2017; Arcavi et al. 2017; Coulter et al. 2017; Lipunov et al. 2017; Yang et al. 2017; Tanvir et al. Levan 2017) with a luminosity, thermal spectrum, and rapid temporal decay consistent with those predicted for "kilonova" (KN) emission, powered by the radioactive decay of heavy elements synthesized in the merger ejecta (Li & Paczyński 1998; Metzger et al. 2010). The presence of both early-time visual ("blue") emission (Metzger et al. 2010), which transitioned to near-infrared ("red") emission (Barnes & Kasen 2013; Tanaka & Hotokezaka 2013) at late times, requires at least two distinct ejecta components consisting, respectively, of light and heavy r-process nuclei (e.g., Chornock et al. 2017; Cowperthwaite et al. 2017; Drout et al. 2017; Kasen et al. 2017; Kasliwal et al. 2017; Nicholl et al. 2017). Rising X-ray (Margutti et al. 2017; Troja et al. 2017) and radio (Alexander et al. 2017; Hallinan et al. 2017) emission was observed roughly two weeks after the merger, consistent with the delayed onset of the synchrotron afterglow of a more powerful relativistic GRB, with an emission that was initially relativistically beamed away from our line of sight (e.g., van Eerten & MacFadyen 2011).

The discovery of GW170817 implies a BNS rate of ${{ \mathcal R }}_{\mathrm{BNS}}={1540}_{-1220}^{+3200}\,{\mathrm{Gpc}}^{-3}\,{\mathrm{yr}}^{-1}$, corresponding to ≈6–120 BNS mergers per year once LIGO/Virgo reach design sensitivity (LIGO Scientific Collaboration & Virgo Collaboration 2017). This relatively high rate bodes well for the prospects of several scientific objectives requiring a large population of GW detections, such as "standard siren" measurements of the cosmic expansion history (Holz & Hughes 2005; Nissanke et al. 2010; LIGO Scientific Collaboration et al. 2017c) or as probes of the equation of state (EOS) of NSs (e.g., Read et al. 2009; Hinderer et al. 2010; Bauswein & Janka 2012).

Uncertainties in the EOS limit our ability to predict key properties of NSs, such as their radii and maximum stable mass (e.g., Özel & Freire 2016). Methods to measure NS radii from GWs include searching for tidal effects on the waveform during the final stages of the BNS inspiral (Damour & Nagar 2010; Hinderer et al. 2010; Damour et al. 2012; Del Pozzo et al. 2013; Read et al. 2013; Favata 2014; Agathos et al. 2015; Chatziioannou et al. 2015; Lackey & Wade 2015) and for quasi-periodic oscillations of the post-merger remnant (e.g., Bauswein & Janka 2012; Bauswein et al. 2012; Clark et al. 2014; Bauswein & Stergioulas 2015; Bauswein et al. 2016). Searches on timescales of tens of ms to $\lesssim 500\,{\rm{s}}$ post-merger revealed no evidence for such quasi-periodic oscillations in GW170817 (LIGO Scientific Collaboration & Virgo Collaboration 2017).

While the radii of NS are controlled by the properties of the EOS at approximately twice the nuclear saturation density, the maximum stable mass (Lattimer & Prakash 2001), ${M}_{\max }$ instead depends on the very high-density EOS (around eight times the saturation density; Özel & Psaltis 2009). Observations of two pulsars with gravitational masses of $1.93\pm 0.07\,{M}_{\odot }$ (Demorest et al. 2010; Özel & Freire 2016) or $2.01\pm 0.04\,{M}_{\odot }$ (Antoniadis et al. 2013) place the best current lower bounds. However, other than the relatively unconstraining limit set by causality, no firm theoretical or observational upper limits exist on ${M}_{\max }$. Indirect, assumption-dependent limits on ${M}_{\max }$ exist from observations of short GRBs (e.g., Lasky et al. 2014; Fryer et al. 2015; Lawrence et al. 2015; Piro et al. 2017) and by modeling the mass distribution of NSs (e.g., Antoniadis et al. 2016; Alsing et al. 2017).

Despite the large uncertainties on ${M}_{\max }$, it remains one of the most important properties affecting the outcome of a BNS merger and its subsequent electromagnetic (EM) signal (Figure 1). If the total binary mass ${M}_{\mathrm{tot}}$ exceeds a critical threshold of ${M}_{\mathrm{th}}\approx {{kM}}_{\max }$, then the merger product undergoes "prompt" dynamical-timescale collapse to a black hole (BH; e.g., Shibata 2005; Shibata & Taniguchi 2006; Baiotti et al. 2008; Hotokezaka et al. 2011), where the proportionality factor k ≈ 1.3−1.6 is greater for smaller values of the NS "compactness", ${C}_{\max }=({{GM}}_{\max }/{c}^{2}{R}_{1.6})$, where ${R}_{1.6}$ is the radius of a 1.6 M_⊙ NS (e.g., Bauswein et al. 2013). For slightly less-massive binaries with ${M}_{\mathrm{tot}}\lesssim {M}_{\mathrm{th}}$, the merger instead produces a hyper-massive neutron star (HMNS), which is supported from collapse by differential rotation (and, potentially, by thermal support). For lower values of ${M}_{\mathrm{tot}}\lesssim 1.2\,{M}_{\max }$, the merger instead produces a supramassive neutron star (SMNS), which remains stable even once its differential rotation is removed, as is expected to occur $\lesssim 10\mbox{--}100$ ms following the merger (Baumgarte et al. 2000; Paschalidis et al. 2012; Kaplan et al. 2014). A SMNS can survive for several seconds, or potentially much longer, until its rigid body angular momentum is removed through comparatively slow processes, such as magnetic spin-down. Finally, for an extremely low binary mass, ${M}_{\mathrm{tot}}\lesssim {M}_{\max }$, the BNS merger produces an indefinitely stable NS remnant (e.g., Bucciantini et al. 2012; Giacomazzo & Perna 2013). Figure 2 shows the baryonic mass thresholds of these possible BNS merger outcomes (prompt collapse, HMNS, SMNS, stable) for an example EOS as vertical dashed lines.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The different types of merger outcomes are predicted to create qualitatively different electromagnetic (EM) signals (e.g., Bauswein et al. 2013; Metzger & Fernández 2014). In this Letter, we combine EM constraints on the type of remnant that formed in GW170817 with GW data on the binary mass in order to constrain the radii and maximum mass of NSs.

This section reviews what constraints can be placed from EM observations on the energy imparted by a long-lived NS into the non-relativistic KN ejecta (Section 2.1) and into the relativistic ejecta of the GRB jet (Section 2.2). Then in Section 2.3 we describe the implications for the type of remnant formed.

2.1. Kilonova (Non-relativistic Ejecta)

2.2. Gamma-Ray Burst (Relativistic Ejecta)

2.3. Constraints on the Merger Remnant in GW170817

The KN emission from GW170817 tightly constrains the type of compact remnant that formed in the merger event (Figure 1). Prompt collapse to a BH (${M}_{\mathrm{tot}}\gtrsim {M}_{\mathrm{th}}$) is disfavored by the quantity of the blue KN ejecta. General relativistic numerical simulations show that mergers with prompt collapses eject only a small quantity $\lesssim {10}^{-4}\mbox{--}{10}^{-3}\,{M}_{\odot }$ of matter from the merger interface (e.g., Hotokezaka et al. 2011), inconsistent with the inferred value ${M}_{\mathrm{ej}}^{\mathrm{blue}}\gtrsim {10}^{-2}\,{M}_{\odot }$ for GW170817. The accretion disk outflows can also contribute; however, the wind ejecta with ${Y}_{e}\gtrsim 0.25$ is only a fraction of the initial torus, which is already small $\lesssim 0.01\mbox{--}0.02\,{M}_{\odot }$ for prompt collapse (Ruffert & Janka 1999; Shibata & Taniguchi 2006; Oechslin et al. 2007b). Furthermore, the predicted velocities of the disk winds $\sim 0.03\mbox{--}0.1c$ (e.g., Fernández & Metzger 2013; Just et al. 2015) are lower than the velocities $\gtrsim 0.2\mbox{--}0.3c$ inferred for the blue KN of GW170817 (e.g., Nicholl et al. 2017).

An HMNS remnant, due to its longer lifetime $\gtrsim 10$ ms, produces a higher quantity of dynamical and disk-wind ejecta. The expansion rate $\sim 0.2\mbox{--}0.3c$ and ejecta mass of ${M}_{\mathrm{ej}}^{\mathrm{blue}}\approx 0.01\mbox{--}0.02\,{M}_{\odot }$ of the blue KN ejecta inferred for GW170817 are consistent with the properties of the high-Y_e shock-heated dynamical ejecta found by BNS merger simulations (e.g., Sekiguchi et al. 2016), provided that the radius of the NS is relatively small, ${R}_{\mathrm{ns}}\lesssim 11$ km (Nicholl et al. 2017). The higher quantity and lower velocity of the red KN emission are also broadly consistent with those expected from the outflows of a relatively massive accretion torus $\approx 0.1\mbox{--}0.2\,{M}_{\odot }$ (e.g., Siegel & Metzger 2017b) following the collapse of a relatively short-lived HMNS.

At the other extreme, a long-lived SMNS or indefinitely stable NS remnant is strongly disfavored by the moderate kinetic energy of the observed KN and GRB afterglow. Even once its differential rotation has been removed, a SMNS possesses an enormous rotational energy, $T\approx {10}^{53}$ erg, which is available to be deposited into the post-merger environment. Not all of this energy is "extractable" insofar as, even just prior to spinning down to the threshold for collapse, the NS remnant is still rotating quite rapidly. Margalit et al. (2015) showed that the collapse of an SMNS is unlikely to produce a centrifugally supported accretion disk outside of the innermost stable circular orbit, in which case all of the mass and angular momentum of the star are trapped in the BH.

The extractable rotational energy of a BNS merger remnant is more precisely defined as

$\begin{eqnarray}&&{\rm{\Delta }}T={T}_{0}-{T}_{\bullet },\end{eqnarray} \tag{ 1 }$

where T₀ is the energy available immediately after differential rotation has been removed, and T_• is the rotational energy at the point of gravitational collapse to a BH. We take T₀ equal to the rotational energy at the mass-shedding limit, a condition that approximates the state of the remnant immediately after differential rotation is removed. However, the constraints obtained hereafter would be similar if we had instead taken T₀ to equal the threshold value $T/| W| \approx 0.14$ for the growth of secular instabilities (e.g., Lai & Shapiro 1995). Figure 2 shows that ${\rm{\Delta }}T$ rises sharply from zero at the HMNS–SMNS boundary to ${\rm{\Delta }}T={T}_{0}\approx {10}^{53}$ erg for stable remnants.

The most likely mechanism by which ${\rm{\Delta }}T$ is removed, enabling the SMNS to collapse, is the extraction of angular momentum via a magnetized outflow or jet. MHD BNS merger simulations find that ultra-strong magnetic fields $\gtrsim {10}^{15}\mbox{--}{10}^{16}$ G are generated in the merger remnant (e.g., Kiuchi et al. 2014). An NS of radius ${R}_{\mathrm{ns}}$, rotation frequency ${\rm{\Omega }}=2\pi /P$, and spin period P loses energy to a magnetic wind at a rate (Spitkovsky 2006)

$\begin{eqnarray}&&{\dot{E}}_{\mathrm{mag}}=\displaystyle \frac{{\mu }^{2}{{\rm{\Omega }}}^{4}}{{c}^{3}}(1+{\sin }^{2}\chi ),\end{eqnarray} \tag{ 2 }$

where ${B}_{{\rm{d}}}$, $\mu ={B}_{{\rm{d}}}{R}_{\mathrm{ns}}^{3}$, and χ is the surface magnetic dipole field strength, dipole moment, and angle between the rotation and dipole axes, respectively. Taking ${R}_{\mathrm{ns}}=12$ km and $\chi =0$, the SMNS's available rotational energy is removed by magnetic torques on a timescale

$\begin{eqnarray}&&{\tau }_{\mathrm{sd}}=\displaystyle \frac{{\rm{\Delta }}T}{{\dot{E}}_{\mathrm{mag}}}\approx 24\,{\rm{s}}\left(\displaystyle \frac{{\rm{\Delta }}T}{{10}^{52}\,\mathrm{erg}}\right){\left(\displaystyle \frac{{B}_{{\rm{d}}}}{{10}^{15}{\rm{G}}}\right)}^{-2}{\left(\displaystyle \frac{P}{0.8\mathrm{ms}}\right)}^{4}.\end{eqnarray} \tag{ 3 }$

If we demand that BH formation occur on a timescale of ${\tau }_{\mathrm{sd}}\lesssim 2\,{\rm{s}}$ following the merger in order to explain the observed gamma-ray emission (Section 2.2), then this requires a SMNS remnant with ${B}_{{\rm{d}}}\gg {10}^{15}$ G or ${\rm{\Delta }}T\ll {10}^{53}$ erg.

An SMNS can in principle also spin down through gravitational-wave emission, as may result from the quadrupolar moment of inertia induced by a strong interior magnetic field that is misaligned with the rotation axis (e.g., Stella et al. 2005; Dall'Osso et al. 2009, 2015). Figure 3 shows that GW spin-down dominates over magnetic spin-down (Equation (2)) only if the interior toroidal magnetic field exceeds the external poloidal one by a factor of $\gtrsim 100$. However, such a strong toroidal to poloidal field configuration would be unstable (Braithwaite 2009; Akgün et al. 2013; gray shaded region in Figure 3) and would furthermore imply a relatively long collapse time of ${\tau }_{\mathrm{sd}}\gtrsim 100\,{\rm{s}}$, potentially incompatible with the gamma-ray burst emission observed on a timescale $\lesssim 2\,{\rm{s}}$ (Section 2.2). It would also produce quasi-periodic GW emission, which is not observed in the GW170817 post-merger signal (albeit with only weakly constraining upper limits; LIGO Scientific Collaboration et al. 2017).

Zoom In Zoom Out Reset image size

In summary, all signs point to GW170817 having produced an HMNS or very short-lived SMNS remnant. If the merger had instead produced a long-lived SMNS, then a large fraction of its available rotational energy $\gtrsim {10}^{52}$ erg should have been deposited into the merger environment, either into a collimated relativistic jet or shared more equitably with the merger ejecta, on a timescale $\sim {\tau }_{\mathrm{sd}}$. Such a large energy input is incompatible with the GRB and KN observations of GW170817.

The masses of the binary components inferred for GW170817, combined with evidence from the KN disfavoring a prompt collapse, places a lower limit on the maximum mass ${M}_{\max }$ of a (slowly rotating) NS. Likewise, upper limits on the rotational energy injected by a long-lived SMNS place an upper limit on ${M}_{\max }$.

In order to translate GW+EM inferences into constraints on the properties of NSs, we use the RNS code (Stergioulas & Friedman 1995) to construct general relativistic rotating hydrostationary NS models for a range of nuclear EOSs. We use the piecewise polytropic approximations to EOSs available in the literature provided in Read et al. 2009 (circles in Figure 4; summarized in Table 1), supplemented by EOSs constructed from the two-parameter piecewise polytropic parameterization of Margalit et al. (2015; triangles in Figure 4). Our simplified parameterization is limited in its ability to model micro-physically motivated EOS with high accuracy, yet it allows us to efficiently survey the EOS parameter space, and we leave it to future work to extend this first analysis with more flexible EOS parameterizations (e.g., Raithel et al. 2016). For each EOS, the uncertainty range of the GW170817-measured binary gravitational mass (LIGO Scientific Collaboration & Virgo Collaboration 2017) translates into a corresponding uncertainty range of baryonic mass, defined by the probability distribution

$\begin{eqnarray}&&P({M}_{\mathrm{rem}}^{{\rm{b}}}| { \mathcal O },\mathrm{EoS})\\ &&\quad =\displaystyle \int {{dM}}_{1}^{{\rm{b}}}\displaystyle \int {{dM}}_{2}^{{\rm{b}}}\,\,\delta ({M}_{1}^{{\rm{b}}}+{M}_{2}^{{\rm{b}}}-{M}_{\mathrm{ej}}-{M}_{\mathrm{rem}}^{{\rm{b}}})\\ &&\times \,P({g}_{\mathrm{EoS}}({M}_{1}^{{\rm{b}}}),{g}_{\mathrm{EoS}}({M}_{2}^{{\rm{b}}})| { \mathcal O })| {g}_{\mathrm{EoS}}^{\backprime}({M}_{1}^{{\rm{b}}})| | {g}_{\mathrm{EoS}}^{\backprime}({M}_{2}^{{\rm{b}}})| .\end{eqnarray} \tag{ 4 }$

Here $P({M}_{1}^{{\rm{g}}},{M}_{2}^{{\rm{g}}}| { \mathcal O })$ is the posterior joint probability distribution function of NS gravitational masses inferred from the BNS waveform ${ \mathcal O }$ (LIGO Scientific Collaboration & Virgo Collaboration 2017), ${M}_{\mathrm{ej}}=2\times {10}^{-2}\,{M}_{\odot }$ is a conservative lower limit for the mass loss from the system as inferred from the KN ejecta, and the EOS enters in converting between gravitational and baryonic masses, ${M}^{{\rm{g}}}={g}_{\mathrm{EoS}}({M}^{{\rm{b}}})$. We approximate the posterior by changing variables to the chirp mass ${{ \mathcal M }}_{{\rm{c}}}$ and mass-ratio $q={M}_{1}^{{\rm{g}}}/{M}_{2}^{{\rm{g}}}$,

$\begin{eqnarray}&&P({M}_{1}^{{\rm{g}}},{M}_{2}^{{\rm{g}}}| { \mathcal O })=P(q,{{ \mathcal M }}_{{\rm{c}}}){{ \mathcal M }}_{{\rm{c}}}^{-1}{q}^{6/5}{(1+q)}^{-2/5},\end{eqnarray} \tag{ 5 }$

assuming independent asymmetric Gaussian distributions for both $P({{ \mathcal M }}_{{\rm{c}}})$ and P(q), consistent with the median and 90% quoted confidence levels on ${{ \mathcal M }}_{{\rm{c}}}$, ${M}_{1}^{{\rm{g}}}$, ${M}_{2}^{{\rm{g}}}$, and ${M}_{\mathrm{tot}}^{{\rm{g}}}$. Specifically, we assume $P(q,{{ \mathcal M }}_{{\rm{c}}})$ $\propto \,\exp [-{({{ \mathcal M }}_{{\rm{c}}}-{\mu }_{{ \mathcal M }},\pm )}^{2}/2{\sigma }_{{ \mathcal M }}^{2}$ $-\,{(q-{\mu }_{q})}^{2}/2{\sigma }_{q}^{2}]$ for $q\leqslant 1$ and P = 0 otherwise, with ${\mu }_{q}=1,{\sigma }_{q}\simeq 0.164,{\mu }_{{ \mathcal M }}\simeq 1.188\,{M}_{\odot }$ and where ${\sigma }_{{ \mathcal M },\pm }\,\simeq 2.63\times {10}^{-3}\,{M}_{\odot }$ ($2.07\times {10}^{-3}\,{M}_{\odot }$) for ${{ \mathcal M }}_{{\rm{c}}}\geqslant {\mu }_{{ \mathcal M }}$ (${{ \mathcal M }}_{{\rm{c}}}\,\lt {\mu }_{{ \mathcal M }}$), respectively.

Zoom In Zoom Out Reset image size

Table 1.  EOS Properties and Consistency with EM Observations

	${M}_{\max }^{{\rm{g}}}$	${R}_{1.3}$	${M}_{\mathrm{smns}}^{{\rm{g}}}$	${\rm{\Delta }}{T}_{\max }$	Consistency
EOS	$({M}_{\odot })$	(km)	$({M}_{\odot })$	(1053erg)	(%)
MS1	2.77	14.9	3.31	1.8	0.0
MPA1	2.45	12.4	2.97	1.8	0.0
APR3	2.37	12.0	2.84	1.7	0.2
ENG	2.24	12.0	2.67	1.4	5.2
WFF2	2.20	11.1	2.63	1.6	10.2
APR4	2.19	11.3	2.61	1.5	18.4
SLy	2.05	11.8	2.43	1.2	100.0
H4	2.02	14.0	2.38	0.8	100.0
ALF2	1.98	12.7	2.41	0.9	100.0
GNH3a	1.96	14.3	2.29	0.7	100.0
ALF4a	1.93	11.5	2.35	1.0	99.8
BBB2a	1.92	11.2	2.27	1.1	99.4
MS2a	1.80	14.3	2.10	0.6	99.9

Note. All EOSs are approximated as piecewise broken polytropes (Read et al. 2009).

^aRuled out by $2.01\pm 0.04\,{M}_{\odot }$ mass of PSR J0348+0432 (Antoniadis et al. 2013).

Download table as:  ASCIITypeset image

For each EOS, we then compare the inferred remnant mass to the "allowed" range between the maximum mass to avoid prompt collapse (using the relation ${M}_{\mathrm{th}}({R}_{1.6},{M}_{\max })$ of Bauswein et al. 2013), to the minimum baryonic mass, which results in an SMNS with an extractable energy ${\rm{\Delta }}T$ (Equation (1)) less than the upper limits on the kinetic energy of the KN and GRB emission ${E}_{\mathrm{EM}}={E}_{\mathrm{KN}}+{E}_{\mathrm{GRB}}\lesssim {10}^{51}\,\mathrm{erg}$. Integrating the probability distribution of the remnant mass within this allowed range yields the "consistency" of the given EOS with the GW170817 observations,

$\begin{eqnarray}&&\mathrm{Consistency}={\int }_{S}P({M}_{\mathrm{rem}}^{{\rm{b}}}| { \mathcal O },\mathrm{EoS})\,{{dM}}_{\mathrm{rem}}^{{\rm{b}}},\end{eqnarray} \tag{ 6 }$

where S is the domain in which both ${\rm{\Delta }}T({M}_{\mathrm{rem}}^{{\rm{b}}})\leqslant {E}_{\mathrm{EM}}$ and ${M}_{\mathrm{rem}}^{{\rm{b}}}\leqslant {M}_{\mathrm{th}}$.

One example of this analysis is illustrated in Figure 2. Clearly, ${E}_{\mathrm{EM}}$ is so much smaller than ${\rm{\Delta }}{T}_{\max }$ that the extractable energy curve intersects ${E}_{\mathrm{EM}}$ at the very precipice of the SMNS–HMNS transition. We also find for all of the EOSs that we have examined that ${M}_{\mathrm{smns}}^{{\rm{b}}}\approx 1.18{M}_{\max }^{{\rm{b}}}$ largely irrespective of compactness, consistent with previous findings (e.g., Lasota et al. 1996). These two facts allow for the formulation of an approximate analytic criterion on the maximal non-rotating NS mass consistent with GW170817,

$\begin{eqnarray}&&{M}_{\max }^{{\rm{b}}}\lesssim {M}_{\mathrm{rem}}^{{\rm{b}}}/\xi ,\end{eqnarray} \tag{ 7 }$

where $\xi \simeq 1.16\mbox{--}1.21$ and the EOS is only necessary in translating baryonic to gravitational masses.

In addition to several key properties of each EOS, Table 1 provides the probability that each EOS is consistent with constraints from GW170817. For instance, the very hard MS1, MPA1, and ENG EOS are disfavored, with consistencies of 0.0%, 0.0% and 5.2%, respectively. However, the softer EOSs with ${M}_{\max }^{{\rm{g}}}\lesssim 2.1\mbox{--}2.2\,{M}_{\odot }$ show much higher consistencies. Figure 4 shows where each of our EOSs lie in this ${M}_{\max }\mbox{--}{R}_{1.3}$ parameter space, with the strength of the symbol representing the probability of its consistency with GW170817. Additionally shown are consistency values for polytropic EOSs of the form $P\propto {\rho }^{1+1/n}$ with indices $n=0.5,0.6,0.7$. These define diagonal curves in the ${M}_{\max }^{{\rm{g}}}\mbox{--}{R}_{1.3}$ plane parameterized by the pressure normalization of the polytrope. Regions of large compactness are ruled out by the requirement of causality, ${R}_{\max }\geqslant 2.82{{GM}}_{\max }^{{\rm{g}}}/{c}^{2}$ (dark shaded region; Koranda et al. 1997), which is conservative as ${R}_{1.3}$ is generally larger than the radius of a maximum mass NS, ${R}_{\max }$. A tighter estimate of ${R}_{1.3}\geqslant 3.1{{GM}}_{\max }^{{\rm{g}}}/{c}^{2}$ is therefore also shown (light shaded region). The background gray curve shows the cumulative probability distribution function that the maximum mass ${M}_{\max }$ is less than a given value. This was calculated by marginalizing over the ${R}_{1.3}$ axis and treating the consistency values as points in a probability distribution function. We weight the EOS with ${M}_{\max }^{{\rm{g}}}$ below $2.01\,{M}_{\odot }$ by a Gaussian prior accounting for consistency with the maximum measured pulsar mass of $2.01\pm 0.04\,{M}_{\odot }$ (Antoniadis et al. 2013). Thus, we find ${M}_{\max }^{{\rm{g}}}\lesssim 2.17\,{M}_{\odot }$ at 90% confidence.

Several works have explored the potentially exotic EM signals of BNS mergers in cases when a long-lived SMNS or stable neutron star remnant is formed (Metzger et al. 2008b; Bucciantini et al. 2012; Yu et al. 2013; Metzger & Bower 2014; Metzger & Piro 2014; Gao et al. 2016; Siegel & Ciolfi 2016a, 2016b). However, one of the biggest lessons from GW170817 was the well-behaved nature of its EM emission, under the simplest case of a relatively short-lived HMNS remnant (Shibata & Taniguchi 2006), with an off-axis afterglow (van Eerten & MacFadyen 2011) and KN emission (Metzger et al. 2010) closely resembling "vanilla" theoretical predictions.

Here we have made explicit the argument that the BNS merger GW170817 formed a HMNS. In combination with the GW-measured binary mass, this inferred outcome places upper and lower limits on the maximum NS mass. The lower limit on ${M}_{\max }$ is not constraining compared to those from well-measured pulsar masses, though tighter constraints would be possible by the future detection of similar fast-expanding blue KN ejecta (indicating HMNS formation) from a future BNS merger with a similar observing inclination but a higher binary mass than GW170817. The lack of a luminous blue KN following the short GRB050509b (Bloom et al. 2006; Metzger et al. 2010; Fong et al. 2017) may implicate a high-mass binary and prompt collapse for this event.¹

On the other hand, our upper limits on ${M}_{\max }\lesssim 2.17\,{M}_{\odot }$ (90% confidence limit; Figure 4) are more constraining than the previous weak upper limits from causality, and less model-dependent than other methods (e.g., Lasky et al. 2014; Fryer et al. 2015; Lawrence et al. 2015; Gao et al. 2016; Alsing et al. 2017). A low value of ${M}_{\max }$ has also been suggested based on Galactic NS radius measurements (Özel et al. 2016; Özel & Freire 2016) and would be consistent with the relatively small NS radius $\lesssim 11$ km inferred from modeling the blue KN (Cowperthwaite et al. 2017; Nicholl et al. 2017). Furthermore, the lack of measurable tidal effects in the inspiral of GW170817 similarly imply a small NS radius and thus a low ${M}_{\max }$ (LIGO Scientific Collaboration & Virgo Collaboration 2017). Upper limits on ${M}_{\max }$ will be improved by the future discovery of EM emission from a merger with a lower total mass than GW170817. Conversely, the detection of a substantially brighter afterglow or faster evolving KN emission could instead point to the formation of a long-lived SMNS or stable remnant. The NS masses measured for GW170817 are broadly consistent with being drawn from Galactic NS population, which is well-fit by a Gaussian of mean $\mu =1.32\,{M}_{\odot }$ and standard deviation $\sigma =0.11\,{M}_{\odot }$ (Kiziltan et al. 2013); this hints that the HMNS formation inferred in GW170817 is likely a common—if not the most frequent—outcome of a BNS merger.

A simple analytic estimate of our result can be obtained from Equation (7), using the approximation ${M}_{{\rm{b}}}={M}_{{\rm{g}}}+0.075{M}_{{\rm{g}}}^{2}$ for the relation between baryonic and gravitational masses (Timmes et al. 1996). From this relation, the total baryonic binary mass is constrained, ${M}_{\mathrm{rem}}^{{\rm{b}}}\lesssim {M}_{\mathrm{tot}}^{{\rm{b}}}\lesssim 3.06\,{M}_{\odot }$. A typical value of $\xi \approx 1.18$ then implies that

$\begin{eqnarray}&&{M}_{\max }^{{\rm{g}}}\lesssim \displaystyle \frac{\sqrt{1+0.3{M}_{\mathrm{rem}}^{{\rm{b}}}/\xi }-1}{0.15}\lesssim 2.2\,{M}_{\odot },\end{eqnarray} \tag{ 8 }$

in agreement with our more elaborately calculated result. We stress that the calculation above is intended only as an approximate analytic estimate, and that neither do we use the Timmes et al. (1996) relation nor do we assume a universal value for ξ in our complete analysis (Section 3).

Our approach differs in several respects from similar works constraining ${M}_{\max }$ (e.g., Lawrence et al. 2015; Fryer et al. 2015). These works generally assume (a) that the creation of a GRB implies that a BH formed, and (b) that BH formation necessarily implies either prompt collapse or an HMNS post-merger remnant. The central engine and emission mechanisms of GRBs are still widely debated in the literature, and the validity of the "GRB = BH" assumption remains unclear. Baryon-pollution by neutrino-driven winds launched off a long-lived NS remnant may hinder ultra-relativistic jets (Murguia-Berthier et al. 2014); however, the Lorentz factors of short GRB jets and GRB 170817A in particular are poorly constrained, and there remains room for the possibility that short GRBs may be powered by strongly magnetized rapidly rotating NSs. NS GRB engines have also been suggested on grounds of the "extended" X-ray emission observed after some short GRBs (e.g., Metzger et al. 2008b; Rowlinson et al. 2013), emission that is difficult to interpret within the BH engine model. Furthermore, the peculiar properties of GRB 170817A, accompanying GW170817, although broadly consistent with a normal GRB viewed off-axis (e.g., Margutti et al. 2017), may also point at a difference between this event and cosmological short GRBs (e.g., Gottlieb et al. 2017), necessitating further caution in the GRB modeling and interpretation. Even if a BH did form as a prerequisite to the GRB in this event (i.e., within ∼2 s post-merger), there is nothing a priori preventing this from occurring through the spin-down-induced collapse of an SMNS merger remnant, negating assumption (b) above. Here we have circumvented both assumptions and GRB engine modeling altogether by instead relying on a simple energetic consideration—that an SMNS merger remnant would inevitably release an enormous amount of rotational energy into the surrounding KN ejecta and circumstellar medium. Therefore, only merger remnants with an extractable rotational energy $\lesssim {10}^{51}\,\mathrm{erg}$ are consistent with the energetics inferred from EM observations of GW170817.

Several uncertainties affect our conclusions. Our upper limits on ${M}_{\max }$ implicitly assume that this is the most important parameter controlling the HMNS–SMNS boundary, and that the suite of EOSs that we have taken are sufficiently "representative" in the requisite sense. We cannot obviously exclude the possibility that an alternative EOS could be found with a large ${M}_{\max }$ that would still be consistent with our observational constraints.

Another uncertainty affecting our conclusions is the possibility that in counting the KN and GRB components of the ejecta, we are somehow "missing" substantial additional energy imparted by a putative SMNS remnant to the environment; however, any such hidden ejecta should be at least mildly relativistic and thus tightly constrained by radio synchrotron emission on timescales of months to years following the merger (Metzger & Bower 2014). Yet another uncertainty is the possibility that a SMNS did form, but most of its rotational energy was lost to GW radiation instead of being transferred to the merger ejecta. Though unlikely, we found this would only be possible for remnant lifetimes of $\sim 100\,{\rm{s}}$ (Figure 3). Searches in the GW170817 waveform have revealed no evidence for such signals, although the detectors' decreasing sensitivity at high frequencies currently limits these constraints (LIGO Scientific Collaboration et al. 2017).

Finally, our analysis neglects the effects of thermal pressure on the stability of the SMNS (Kaplan et al. 2014), which can be important on timescales of hundreds of milliseconds to seconds post-merger (depending also on the effects of neutrino-driven convection; Roberts et al. 2012). Thermal pressure in the outer layers of the star generally acts to reduce the maximum mass of the SMNS remnant by up to ≲8% (mainly by reducing the angular velocity at the mass-shedding limit), which would act to weaken our constraints on ${M}_{\max }$. Future numerical work exploring the transition from the HMNS to SMNS phase, which includes the effects of neutrino cooling and convection self-consistently, is required to better understand how this would quantitatively affect our conclusions.

B.M. and B.D.M. are supported in part by NASA through the ATP program, grant numbers NNX16AB30G and NNX17AK43G. The statistical analysis in this work would have greatly benefited from open access to the best-fit posterior parameter distributions obtained by LIGO Scientific Collaboration & Virgo Collaboration (2017).