GW170817: Joint Constraint on the Neutron Star Equation of State from Multimessenger Observations

The properties of matter at supranuclear densities determining the internal structure and mass–radius relation of neutron stars (NSs), are poorly known at the moment (Ozel & Freire 2016). Presently, the strongest constraint comes from the fact that the maximum mass for NSs must be larger than about $2\,{M}_{\odot }$ (Antoniadis et al. 2013). Gravitational wave (GW) observations of coalescing binary NSs have long been considered as a promising avenue to constrain the equation of state (EOS) of dense matter. The tidal polarizability of the NSs is encoded in the phase evolution of the gravitational wave (GW) signal during the inspiral (Flanagan & Hinderer 2008; Damour & Nagar 2010; Hinderer et al. 2010, 2016; Damour et al. 2012; Del Pozzo et al. 2013; Read et al. 2013; Favata 2014; Wade et al. 2014; Bernuzzi et al. 2015b; Lackey & Wade 2015; Hotokezaka et al. 2016; Lackey et al. 2017; Dietrich et al. 2017a; Kiuchi et al. 2017). The post-merger signal, if detected, could also place strong constraints on the physics of high-density matter (Bauswein & Janka 2012; Takami et al. 2014; Bernuzzi et al. 2015a; Chatziioannou et al. 2017; Radice et al. 2017a; Yang et al. 2017).

On 2017 August 17, gravitational waves (GWs) from a pair of merging NSs were observed, for the first time, by the LIGO-Virgo detector network (Abbott et al. 2017b): GW170817. Less than 2 seconds after the end of the GW signal, a short γ-ray burst was detected by the Fermi and INTEGRAL satellites in a coincident sky position (Abbott et al. 2017a). In the following hours and days, the same source, now named AT2017gfo, was detected in the X-ray, UV, optical, infrared, and radio bands (Abbott et al. 2017c; Arcavi et al. 2017; Chornock et al. 2017; Coulter et al. 2017; Cowperthwaite et al. 2017; Drout et al. 2017; Evans et al. 2017; Hallinan et al. 2017; Kasliwal et al. 2017; Murguia-Berthier et al. 2017; Nicholl et al. 2017; Smartt et al. 2017; Soares-Santos et al. 2017; Tanaka et al. 2017; Tanvir et al. 2017; Troja et al. 2017).

The preliminary analysis of GW170817 presented in Abbott et al. (2017b) already provided a first constraint on the amplitude of tidal effects during the binary inspiral, disfavoring EOSs with large NS radii. Margalit & Metzger (2017) argued that the merger remnant might not have formed a long-lived remnant because of the relatively low energy of the ejecta inferred from optical and infrared data. Under this assumption, Margalit & Metzger (2017), and subsequently Shibata et al. (2017), Rezzolla et al. (2017), and Ruiz et al. (2017), placed upper bounds on the maximum mass supported by the NS EOS. Bauswein et al. (2017) pointed out that a prompt black hole (BH) formation is also unlikely, because this would have suppressed the ejection of matter and the subsequent emissions in the optical/infrared. Bauswein et al. (2017) combined this observation with empirical relations between NS radii and the threshold mass for prompt collapse, which was previously found by means of simulations with an approximate treatment of general relativity (GR; Bauswein et al. 2013a), to tentatively rule out EOSs predicting very small NS radii.

In this Letter, we propose and apply to GW170817 a new approach that combines optical/infrared and GW observations, by means of new numerical-relativity results, to derive strong joint constraints on the tidal deformability of NSs.

The GW data tightly constrained the 90% credible interval for the chirp mass of the binary, ${{ \mathcal M }}_{\mathrm{chirp}}\,={({M}_{A}{M}_{B})}^{3/5}{({M}_{A}+{M}_{B})}^{-1/5}$, M_A and M_B being the NS masses, to be ${1.188}_{-0.002}^{+0.004}\,{M}_{\odot }$ (Abbott et al. 2017b). With the same confidence, the binary mass ratio $q={M}_{B}/{M}_{A}$ is constrained to be 0.7–1.0 if the dimensionless NS spins are less than 0.05 (Abbott et al. 2017b). If the priors on the NS spins are relaxed, q becomes only constrained to be within 0.4–1.0. Note, however, that large spins are not expected on the basis of the observed galactic NS binary population (Abbott et al. 2017b). Moreover, $q\lt 0.7$ for this event would imply an implausible mass for the secondary NS, smaller than $1.15\,{M}_{\odot }$, in tension with core-collapse supernova theory (e.g., Radice et al. 2017b). Finally, we remark that the GW data already places strong limits on the component of the NS spin aligned with the orbital angular momentum (Abbott et al. 2017b).

LIGO and Virgo observations also constrain tidal effects in the inspiral by placing an upper bound on the dimensionless quantity (Flanagan & Hinderer 2008; Favata 2014)

$$\begin{eqnarray}&&\tilde{{\rm{\Lambda }}}=\displaystyle \frac{16}{13}\left[\displaystyle \frac{({M}_{A}+12{M}_{B}){M}_{A}^{4}{\tilde{{\rm{\Lambda }}}}_{A}}{{({M}_{A}+{M}_{B})}^{5}}+(A\leftrightarrow B)\right],\end{eqnarray} \tag{ 1 }$$

which is inferred to be smaller than 800 at the 90% confidence level (Abbott et al. 2017b). In the previous equation

$$\begin{eqnarray}&&{\tilde{{\rm{\Lambda }}}}_{i}=\displaystyle \frac{2}{3}{k}_{2}^{(i)}{\left[\left(\displaystyle \frac{{c}^{2}}{G}\right)\left(\displaystyle \frac{{R}_{i}}{{M}_{i}}\right)\right]}^{5},\quad i=A,B\end{eqnarray} \tag{ 2 }$$

are the dimensionless quadrupolar tidal parameters (or tidal polarizability coefficients), where ${k}_{2}^{(i)}$ are the quadrupolar Love numbers for each star. The fate of the merger remnant is not known. The post-merger high-frequency GWs were too weak to be detected, so information on the remnant is not available from GW observations (Abbott et al. 2017d).

The optical and infrared electromagnetic (EM) data is well explained by the radioactive decay of $\sim 0.05\,{M}_{\odot }$ of material (Chornock et al. 2017; Cowperthwaite et al. 2017; Drout et al. 2017; Nicholl et al. 2017; Perego et al. 2017; Rosswog et al. 2017; Tanaka et al. 2017; Tanvir et al. 2017; Villar et al. 2017). UV/optical light-curve modeling of the early emissions, hours to days after merger, points to the presence of a relatively fast, $v\simeq 0.3\,c$, $M\simeq 0.02\,{M}_{\odot }$, component of the outflow (Cowperthwaite et al. 2017; Drout et al. 2017; Nicholl et al. 2017; Perego et al. 2017; Villar et al. 2017). The modeling of the later optical/infrared data points to the presence of at least another component of the outflow with $v\simeq 0.1\,c$ and $M\simeq 0.04\ {M}_{\odot }$ (Chornock et al. 2017; Cowperthwaite et al. 2017; Drout et al. 2017; Perego et al. 2017; Villar et al. 2017). The inferred effective opacities for these two (or more) outflow components suggest that they had different compositions and, possibly, different origins.

GR simulations indicate that only up to $\sim 0.01\,{M}_{\odot }$ of material can be unbound dynamically during the merger itself (Bauswein et al. 2013b; Hotokezaka et al. 2013; Lehner et al. 2016; Radice et al. 2016; Sekiguchi et al. 2016; Bovard et al. 2017; Dietrich et al. 2017b), although larger ejecta masses can be reached for small mass ratios $q\lesssim 0.6$ (Dietrich et al. 2017c). The largest ejecta masses are obtained for soft EOSs. In these cases, the outflows are fast, $v\simeq (0.2\mbox{--}0.4)c$, shock heated, and re-processed by neutrinos (Sekiguchi et al. 2015; Foucart et al. 2016; Radice et al. 2016). Consequently, the dynamic ejecta can potentially explain the UV/optical emissions in the first hours to days. The inferred properties for the outflow component powering the optical/infrared emission on a days to weeks timescale are more easily explained by neutrino, viscous, or magnetically driven outflows from the merger remnant (Metzger et al. 2008, 2009; Dessart et al. 2009; Fernández & Metzger 2013; Metzger & Fernández 2014; Perego et al. 2014; Siegel et al. 2014; Just et al. 2015; Wu et al. 2016; Lippuner et al. 2017; Siegel & Metzger 2017). Detailed modeling suggests that a disk mass of at least $0.08\,{M}_{\odot }$ is required to explain AT2017gfo (Perego et al. 2017).

We perform 29 merger simulations using the GR hydrodynamics code WhiskyTHC (Radice & Rezzolla 2012; Radice et al. 2014a, 2014b). We consider both equal and unequal mass configurations, and we adopt four temperature and composition dependent nuclear EOSs spanning the range of the nuclear uncertainties: the DD2 EOS (Hempel & Schaffner-Bielich 2010; Typel et al. 2010), the BHBΛϕ EOS (Banik et al. 2014), the LS220 EOS (Lattimer & Swesty 1991), and the SFHo EOS (Steiner et al. 2013). This is the largest data set of simulations performed in full-GR and with realistic microphysics to date. Neutrino cooling and Y_e evolution are treated as discussed in Radice et al. (2016). The computational setup is the same as in Radice et al. (2017a). The resolution of the grid regions covering the NSs and the merger remnant is $\simeq 185\,{\rm{m}}$. We verify the robustness of our results and estimate the numerical uncertainties by performing six additional simulations at 25% higher resolution. We conservatively estimate finite-resolution error on the disk and dynamic ejecta masses to be

$$\begin{eqnarray}&&{\rm{\Delta }}{M}_{\mathrm{disk},\mathrm{ej}}=0.5\,{M}_{\mathrm{disk},\mathrm{ej}}+{\epsilon }_{\mathrm{disk},\mathrm{ej}},\end{eqnarray} \tag{ 3 }$$

where ${\epsilon }_{\mathrm{disk}}=5\times {10}^{-4}{M}_{\odot }$ and ${\epsilon }_{\mathrm{ej}}=5\times {10}^{-5}{M}_{\odot }$. A more detailed account of these simulations will be given elsewhere (D. Radice et al. 2017, in preparation). A summary of the simulations is given in Table 1.

Table 1.  Gravitational and Baryonic Masses, Compactnesses, Tidal Deformability Parameters, BH Formation Time, and Disk and Ejecta Masses

EOS	MAa	MBa	${M}_{A}^{* }$ b	${M}_{B}^{* }$ b	CAc	CBc	${\tilde{{\rm{\Lambda }}}}_{A}$ d	${\tilde{{\rm{\Lambda }}}}_{B}$ d	$\tilde{{\rm{\Lambda }}}$ e	${M}_{\mathrm{disk}}$ f	${M}_{\mathrm{ej}}$ g	${t}_{\mathrm{BH}}$ h	${t}_{\mathrm{end}}$ i
	$({M}_{\odot })$	$({M}_{\odot })$	$({M}_{\odot })$	$({M}_{\odot })$						$({10}^{-2}\,{M}_{\odot })$		$(\mathrm{ms})$	$(\mathrm{ms})$
BHBΛϕ	1.365	1.25	1.491	1.352	0.153	0.140	805	1310	1028	18.73	0.06	...	23.98
BHBΛϕ	1.35	1.35	1.473	1.473	0.151	0.151	857	857	857	14.45	0.07	...	21.26
BHBΛϕ	1.4	1.2	1.533	1.297	0.157	0.135	697	1630	1068	20.74	0.11	...	23.74
BHBΛϕ	1.4	1.4	1.533	1.533	0.157	0.157	697	697	697	7.05	0.09	11.96	16.39
BHBΛϕ	1.44	1.39	1.580	1.520	0.161	0.155	591	726	655	8.28	0.06	10.39	15.77
BHBΛϕ	1.5	1.5	1.657	1.657	0.168	0.168	462	462	462	1.93	0.05	2.27	11.78
BHBΛϕ	1.6	1.6	1.778	1.778	0.179	0.179	306	306	306	0.09	0.00	0.99	10.67
DD2	1.365	1.25	1.491	1.352	0.153	0.140	807	1309	1028	20.83	0.04	...	24.24
DD2	1.35	1.35	1.473	1.473	0.151	0.151	858	858	858	15.69	0.03	...	24.41
DD2	1.4	1.2	1.533	1.297	0.157	0.135	699	1630	1070	19.26	0.09	...	23.59
DD2	1.4	1.4	1.533	1.533	0.157	0.157	699	699	699	12.36	0.04	...	24.52
DD2	1.44	1.39	1.580	1.520	0.161	0.155	595	728	658	14.40	0.05	...	23.52
DD2	1.5	1.5	1.657	1.657	0.167	0.167	469	469	469	16.70	0.07	...	23.12
DD2	1.6	1.6	1.778	1.778	0.178	0.178	317	317	317	1.96	0.12	2.28	12.08
LS220	1.2	1.2	1.296	1.296	0.139	0.139	1439	1439	1439	17.43	0.14	...	23.22
LS220	1.365	1.25	1.491	1.355	0.159	0.145	636	1119	848	16.86	0.11	...	26.71
LS220	1.35	1.35	1.473	1.473	0.157	0.157	684	684	684	7.25	0.06	20.34	23.84
LS220	1.4	1.2	1.535	1.296	0.163	0.139	536	1439	893	22.82	0.19	...	23.52
LS220	1.4	1.4	1.535	1.535	0.163	0.163	536	536	536	4.58	0.14	9.93	26.95
LS220	1.44	1.39	1.581	1.520	0.168	0.162	442	563	499	3.91	0.19	7.22	14.83
LS220	1.45	1.45	1.596	1.596	0.169	0.169	421	421	421	2.05	0.16	2.26	11.83
LS220	1.6	1.6	1.790	1.790	0.189	0.189	202	202	202	0.07	0.03	0.63	10.42
LS220	1.71	1.71	1.928	1.928	0.205	0.205	116	116	116	0.06	0.03	0.49	9.94
SFHo	1.365	1.25	1.504	1.364	0.169	0.155	393	680	520	8.81	0.15	...	26.41
SFHo	1.35	1.35	1.486	1.486	0.167	0.167	422	422	422	6.23	0.35	11.96	22.88
SFHo	1.4	1.2	1.547	1.303	0.174	0.148	334	868	546	11.73	0.12	...	24.31
SFHo	1.4	1.4	1.547	1.547	0.174	0.174	334	334	334	0.01	0.04	1.07	13.91
SFHo	1.44	1.39	1.598	1.535	0.179	0.173	277	350	312	0.09	0.04	0.87	7.06
SFHo	1.46	1.46	1.623	1.623	0.182	0.182	252	252	252	0.02	0.00	0.70	9.51

Notes. Disk and ejecta masses are given at the final simulation time.

^aNS gravitational mass. ^bNS baryonic mass. ^cNS compactness, ${GM}/R\,{c}^{2}$. ^dDimensionless quadrupolar tidal parameters, Equation (2). ^eDimensionless tidal parameter, Equation (1). ^fGravitationally bound material with $\rho \leqslant {10}^{13}\ {\rm{g}}\ {\mathrm{cm}}^{-3}$ outside of the apparent horizon. ^gDynamic ejecta mass, computed as from the flux of unbound matter through the coordinate-sphere $r=443\ \mathrm{km}$. ^hBH formation time, in milliseconds after merger. ⁱFinal simulation time, in milliseconds after merger.

Download table as:  ASCIITypeset image

We compute the mass of the dynamic ejecta and of the remnant accretion disk for each model. Our results are shown in Table 1 and Figure 1. The typical dynamic ejecta masses in our simulations are of the order of $\sim {10}^{-3}\,{M}_{\odot }$, in good qualitative agreement with previous numerical-relativity results. We do not find any clear indication of a trend in the dynamic ejecta masses as a function of the binary parameters or EOS. However, we find a clear correlation between the disk masses and the tidal parameter $\tilde{{\rm{\Lambda }}}$. According to our simulations, binaries with $\tilde{{\rm{\Lambda }}}\lesssim 450$ inevitably produce BHs with small $\lesssim {10}^{-2}\,{M}_{\odot }$ accretion disks. These cases are incompatible with the infrared data for AT2017gfo, even under the assumption that all of the matter left outside of the event horizon will be ejected.

Zoom In Zoom Out Reset image size

The reason for this trend is easily understood from the lower panel of Figure 1. The NS dimensionless quadrupolar tidal parameters depend on the negative-fifth power of the NS compactness (${GM}/R\,{c}^{2};$ Equation (2)). Consequently, small values of $\tilde{{\rm{\Lambda }}}$ are associated with binary systems having compact NSs that result in rapid or prompt BH formation. In these cases, the collapse happens on a shorter timescale than the hydrodynamic processes responsible for the formation of the disk. Consequently, only a small amount of mass is left outside of the event horizon at the end of the simulations.

Binaries with larger values of $\tilde{{\rm{\Lambda }}}$ produce more massive disks, up to $\sim 0.2\,{M}_{\odot }$, and longer lived remnants. In these cases, neutrino driven winds and viscous and magnetic processes in the disk are expected to unbind sufficient material to explain the optical and infrared observations for AT2017gfo (Perego et al. 2014; Wu et al. 2016; Siegel & Metzger 2017).

On the basis of our simulations and the current interpretation of the UV/optical/infrared data we can conclude that values of $\tilde{{\rm{\Lambda }}}$ smaller than 400 are tentatively excluded. Together with the LIGO-Virgo constraints on $\tilde{{\rm{\Lambda }}}$ (Abbott et al. 2017b), this result already yields a strong constraint on the EOS.

To illustrate this, we notice that, since the chirp mass of the binary progenitor of GW170817 is well measured, for any given EOS the predicted $\tilde{{\rm{\Lambda }}}$ reduces to a simple function of the mass ratio, that is,

$$\begin{eqnarray}&&\tilde{{\rm{\Lambda }}}=\tilde{{\rm{\Lambda }}}(q,{{ \mathcal M }}_{\mathrm{chirp}}=1.188\,{M}_{\odot };\mathrm{EOS}).\end{eqnarray} \tag{ 4 }$$

We consider a set of 12 EOSs: the four used in the simulations and the other eight from Read et al. (2009). We compute $\tilde{{\rm{\Lambda }}}(q)$ for each and show the resulting curves in Figure 2. There, we also show the upper bound on $\tilde{{\rm{\Lambda }}}$ from the GW observations as well as the newly estimated lower bound from the EM data. On the one hand, stiff EOSs, such as H4 and HB, are already disfavored on the basis of the GW data alone. On the other hand, EOSs as soft as FPS and APR4 are also tentatively excluded on the basis of the EM observations.⁶ Soft EOSs commonly used in simulations, such as SFHo and SLy, lay at the lower boundary of the allowed region, while DD2 and BHBΛϕ are on the upper boundary.

Zoom In Zoom Out Reset image size

Our results show that numerical-relativity simulations are key to exploiting the potential of multimessenger observations. While GW data bounds the tidal deformability of NSs from above, the EM data and our simulations bound it from below. The result is a competitive constraint already after the first detection of a merger event. Our method is general, it can be applied to future observations and used to inform the priors used in the GW data analysis. We anticipate that, with more observations and more precise simulations, the bounds on the tidal deformability of NSs will be further improved.

The physics setting the lower bound on $\tilde{{\rm{\Lambda }}}$ is well understood and under control in our simulations. However, a more extended analysis taking into account the uncertainties in the interpretation of the EM observations and in the simulations is a necessary next step. For example, large components of the NS spins parallel to the orbital plane are not expected, but also not constrained for GW170817. We cannot exclude that, if present, they will affect our results. Moreover, there are indications that small mass ratio binaries $q\lesssim 0.8$ might also form disks with masses up to $\sim 0.1\,{M}_{\odot }$ (Shibata et al. 2017). If confirmed, this would imply that the lower bound on $\tilde{{\rm{\Lambda }}}$ might depend on q. Note that the upper bound on $\tilde{{\rm{\Lambda }}}$ estimated from the GW signal is also likely to have some dependency on q. Consequently, a more precise determination of the exclusion region on $\tilde{{\rm{\Lambda }}}$ will necessarily require a full Bayesian analysis of the GW data using $\tilde{{\rm{\Lambda }}}$ priors informed by numerical-relativity results.

We plan to improve our modeling by means of new simulations exploring the set of binary progenitor parameters compatible with GW170817 and the associated EM counterparts.

It is a pleasure to acknowledge A. Burrows for the many stimulating discussions, and T. Venumadhav for comments on an earlier version of the manuscript. D.R. 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-1523261). D.R. and A.P. acknowledge support from the Institute for Nuclear Theory (17-2b program). S.B. acknowledges support by the EU H2020 under ERC Starting Grant, No. BinGraSp-714626. Computations were performed on the supercomputers Bridges, Comet, and Stampede (NSF XSEDE allocation TG-PHY160025), on NSF/NCSA Blue Waters (NSF PRAC ACI-1440083), Marconi (PRACE proposal 2016153522), and PizDaint/CSCS (ID 667). This manuscript has been assigned LIGO report number LIGO-P1700421 and Virgo report number VIR-0894A-17.