Constraints on the Neutron Star Equation of State from GW170817

TL;DR

The paper reviews how the gravitational wave event GW170817 has been used to constrain the neutron star equation of state, revealing insights into neutron star radii, maximum mass, and the importance of disentangling thermal effects.

Contribution

It provides a comprehensive review of constraints on the neutron star EOS from GW170817 and introduces new radius estimates based on the effective tidal deformability relationship.

Findings

Neutron star radius from GW170817 is between 9.8 and 13.2 km.

Maximum neutron star mass constrained to be less than 2.3 solar masses.

Bayesian inference methods have been used to directly constrain the EOS.

Abstract

The first detection of gravitational waves from a neutron star-neutron star merger, GW170817, has opened up a new avenue for constraining the ultradense-matter equation of state (EOS). The deviation of the observed waveform from a point-particle waveform is a sensitive probe of the EOS controlling the merging neutron stars' structure. In this topical review, I discuss the various constraints that have been made on the EOS in the year following the discovery of GW170817. In particular, I review the surprising relationship that has emerged between the effective tidal deformability of the binary system and the neutron star radius. I also report new results that make use of this relationship, finding that the radius inferred from GW170817 lies between 9.8 and 13.2 km at 90% confidence, with distinct likelihood peaks at 10.8 and 12.3 km. I compare these radii, as well as those inferred in the literature, to X-ray measurements of the neutron star radius. I also summarize the various maximum mass constraints, which point towards a maximum mass < 2.3 M_sun, depending on the fate of the remnant, and which can be used to additionally constrain the high-density EOS. I review the constraints on the EOS that have been performed directly, through Bayesian inference schemes. Finally, I comment on the importance of disentangling thermal effects in future EOS constraints from neutron star mergers.

Equations

$${\mathcal{M}}_c=\frac{(m_1{m}_2{)}^{3/5}}{(m_1+m_2{)}^{1/5}},$$

$$\lambda \equiv -\frac{Q^{(\mathrm{tid})}}{{ϵ}^{(\mathrm{tid})}},$$

$$\Lambda \equiv \frac{\lambda }{M^5}\equiv \frac{2}{3}{k}_2^{(\mathrm{tid})}{C}^{-5},$$

$$\tilde{\Lambda }\equiv \frac{16}{13}\frac{(m_1+12m_2)m_1^4{\Lambda }_1+(m_2+12m_1)m_2^4{\Lambda }_2}{(m_1+m_2{)}^5},$$

$$C=a_0+a_1\ln \Lambda +a_2(\ln \Lambda {)}^2,$$

$$\Lambda (1.4M_{\odot })=2.88\times 10^{-6}{\left(\frac{R(1.4M_{\odot })}{\mathrm{km}}\right)}^{7.5},$$

$${\tilde{\Lambda }}_{qN}={\tilde{\Lambda }}_0\left(1+{\delta }_0(1-q{)}^2\right)+\mathcal{O}\left((1-q{)}^3\right),$$

$${\tilde{\Lambda }}_0=\frac{15-{\pi }^2}{3{\pi }^2}{\xi }^{-5}(1-2\xi {)}^{5/2},$$

$${\delta }_0=\frac{3}{104}(1-2\xi {)}^{-2}\left(-10+94\xi -83{\xi }^2\right),$$

$$\xi =\frac{2^{1/5}G{\mathcal{M}}_c}{R{c}^2}$$

$$P(R)=P(\widetilde{\Lambda})\biggr{\rvert}\frac{\partial\widetilde{\Lambda}}{\partial R}\biggr{\lvert}.$$

$$\Lambda \simeq a{C}^{-6}$$

$$\tilde{\Lambda }\simeq \frac{16a}{13}{\left(\frac{\hat{R}{c}^2}{G{\mathcal{M}}_c}\right)}^{6}f(q),$$

$$f(q)=q^{8/5}(12-11q+12q^2)(1+q{)}^{-26/5}.$$

$$\tilde{\Lambda }=a\prime {\left(\frac{\hat{R}{c}^2}{G{\mathcal{M}}_c}\right)}^6,$$

$$\hat{R}\simeq R_{1.4}\simeq (11.2\pm 0.2)\frac{{\mathcal{M}}_c}{M_{\odot }}{\left(\frac{\tilde{\Lambda }}{800}\right)}^{1/6}\mathrm{km}$$

Full text

1 Introduction

On August 17, 2017, the first detection of gravitational waves from a neutron star-neutron star merger heralded the start of a new era of astrophysics. The event, GW170817, was observed by the two Advanced LIGO and one Advanced Virgo detectors, with a matched-filter signal-to-noise ratio of 32.4. The signal remained in the detectors’ sensitive band for $\sim 100$ s, corresponding to $\sim 3000$ orbital cycles in the initial frequency band considered Abbott2017a . Following the binary coalescence by 1.7 s, the gamma-ray burst (GRB), GRB 170817A, was observed independently by the Fermi Gamma-ray Burst Monitor and the International Gamma-Ray Astrophysics Laboratory, in the same direction as GW170817 Abbott2017c . Roughly 11 hours later, a fading optical transient was discovered by the Swope Supernova Survey (SSS17a/AT2017gfo) Coulter2017 , which was extensively followed up Abbott2017b . Early observations indicated a blue transient that faded over the initial 48 hours, followed by significant reddening over the next $\sim 10$ days. Altogether, these observations indicate that GW170817 was the gravitational signature produced by a neutron star-neutron star merger, which subsequently produced a short-GRB and a kilonova, the latter powered by the radioactive decay of matter in the ejecta (see Ref. Abbott2017b and references therein).

The initial analysis by the LIGO-Virgo collaboration used a matched-filter analysis to determine the parameters of a tidally-corrected post-Newtonian waveform model. The component masses of GW170817 were determined to be in the range $1.17-1.60~{}M_{\odot}$, with a total binary mass of $2.74\begin{subarray}{c}+0.04\\ -0.01\end{subarray}~{}M_{\odot}$. The chirp mass, defined as

[TABLE]

where $m_{1}$ and $m_{2}$ are the component masses, was very tightly constrained, to $\mathcal{M}_{c}=1.188\begin{subarray}{c}+0.004\\ -0.002\end{subarray}~{}M_{\odot}$, for the case of low-spin priors Abbott2018 .111Throughout this paper, I will only report results for the case of low-spin priors, as these are most consistent with the spin distribution inferred for galactic binary neutron stars.

Of particular interest to the neutron star community, was the accompanying first measurement of equation-of-state (EOS) effects on the gravitational wave signal. The waveform produced by the coalescence of two neutron stars should, in principle, deviate from a point-particle waveform, due to tidal effects on the neutron star matter. In one early study, Ref. Read2009a found that these differences could potentially be measured with Advanced LIGO, and that the degree of deviation could be used to differentiate between EOS that differ in the radii they predict for neutron stars by only $\sim 1$ km Read2009a ; Read2013 ; Lackey2015 .

In anticipation of future gravitational wave detections, Ref. Flanagan2008 showed that the EOS-dependent effects during the early phase of the inspiral could be cleanly represented by a single parameter: the tidal Love number, $\lambda$. The tidal Love number is defined as the ratio of the tidally-induced quadrupole deformation, $Q^{(\rm tid)}$, to the tidal potential caused by the binary companion, $\epsilon^{(\rm tid)}$, i.e.,

[TABLE]

or, in its dimensionless form,

[TABLE]

where $M$ is the stellar mass and $C\equiv GM/Rc^{2}$ is the stellar compactness, with $R$ as the radius and $G$ and $c$ the gravitational constant and speed of light, respectively. Here, I follow the convention of Ref. Flanagan2008 and call $k_{2}^{\rm(tid)}$ the tidal apsidal constant. The tidal apsidal constant depends both on the stellar compactness and the particular EOS, and has been constrained to $0.05\lesssim k_{2}^{(\rm tid)}\lesssim 0.15$ for a variety of realistic, hadronic EOS Hinderer2008 ; Hinderer2010 ; Postnikov2010 . Thus, a measurement of $\Lambda$ offers possible insight to the underlying EOS.

The phase of a gravitational-wave signal is determined, to first order, by the related parameter $\widetilde{\Lambda}$, which is the effective tidal deformability of the binary system. The effective tidal deformability is defined as

[TABLE]

where $\Lambda_{1}$ and $\Lambda_{2}$ are the dimensionless tidal deformabilities of the component stars Flanagan2008 ; Favata2014 .

The initial analysis of GW170817 by the LIGO-Virgo Collaboration inferred $\widetilde{\Lambda}\leq 900$ at the 90% confidence level,222The initial analysis in Ref. Abbott2017a quoted a 90% upper limit of $\widetilde{\Lambda}\leq 800$; however, this was an error. The correct upper limit for the original TaylorF2 waveform and low-spin priors is $\widetilde{\Lambda}\leq 900$, as discussed in the caption of Table IV in Ref. Abbott2019 . All of the analyses discussed in this review used the initial upper limit of 800 where relevant. and showed that this disfavors stiff EOS that predict less compact (i.e., more extended) stars Abbott2017a ; Abbott2019 . This measurement prompted a flurry of papers, investigating the implications for the neutron star EOS. Roughly one year after the initial detection, the LIGO-Virgo Collaboration released a re-analysis of the data from GW170817. The re-analysis provided in Ref. Abbott2019 improved on the initial constraints by re-calibrating the Virgo data, extending the range of frequencies included, using a new set of waveform models that go beyond the post-Newtonian approximation, and incorporating source distance measurements from the electromagnetic counterpart. With the new analysis, the chirp mass was revised down slightly to $\mathcal{M}_{c}=1.186\begin{subarray}{c}+0.001\\ -0.001\end{subarray}~{}M_{\odot}$ and more detailed information on the tidal deformability was provided. In particular, they reported a 90% highest posterior density interval of $\widetilde{\Lambda}=300\begin{subarray}{c}+420\\ -230\end{subarray}$.

In this paper, I will review the status of EOS constraints that have been published since the discovery of GW170817 and its initial analysis. I will focus on EOS constraints inferred primarily from the gravitational wave signal, although I note that there is rich information to be learned also from study of the electromagnetic counterpart (see, e.g., Refs. Chornock2017 ; Kasen2017 ; Cote2018 ; Hotokezaka2018 for the implications of GW170817 for r-process nucleosynthesis) as well as from a broader exploration of the gravitational wave data (see, e.g., Refs. Malik2018 ; Zhang2018 ; Carson2019 ; Zhang2019 for constraints on the nuclear symmetry energy or Refs. Montana2018 ; Paschalidis2018 ; Alvarez-Castillo2019 ; Christian2019 for discussion of GW170817 and so-called “twin stars”).

I will start with a brief overview of EOS models and existing constraints. In $\S$3, I will review a series of papers that have found a direct link between $\widetilde{\Lambda}$ and the neutron star radius. In $\S$4, I will review new constraints on the upper limit of the maximum mass of neutron stars, and their implication for the EOS. In $\S$5, I will discuss findings from Bayesian inference schemes that infer EOS parameters directly from the gravitational wave data. Finally, I will comment on the importance of disentangling thermal effects from the cold EOS in future work.

2 Existing equation of state constraints

The EOS of cold, ultra-dense matter remains poorly constrained at high densities, despite decades of study. A wide variety of theoretical models have been proposed, ranging from pure nucleonic models (e.g., Baym1971 ; Friedman1981 ; Akmal1998 ; Douchin2001 ) to models predicting more exotic phases of matter, such as pion condensates (e.g., Pandharipande1975 ), kaon condensates (e.g., Kaplan1986 ), hyperons (e.g., Balberg1997 ), or deconfined quark matter at high densities (e.g., Collins1975 ). More recently, some studies have started to incorporate these quark degrees of freedom using state-of-the-art results from perturbative QCD (e.g., Fraga2014 ). However, despite the diversity of these models, it remains possible that they do not cover the full range of allowed physics. This possibility has motivated the creation of a large number of parametric EOS models. Among these models are the piecewise polytropic or linear representations Read2009 ; Ozel2009 ; Hebeler2010 ; Steiner2010 ; Raithel2016 , spectral expansion methods Lindblom2010 , and constant speed-of-sound parameterizations Alford2013 .

Nuclear experiments are only able to constrain these various models up to densities near the nuclear saturation density, $\rho_{\rm sat}=2.7\times 10^{14}$ g/cm3; extrapolations to higher densities remain difficult Lattimer2012 . Moreover, while the interactions between particles can be written in terms of few-body potentials up to $\sim\rho_{\rm sat}$ Akmal1998 ; Morales2002 ; Gandolfi2012 , the expansion of interactions into few-body terms starts to break down at higher densities, due in part to the overlap of the particle wave-functions Ozel2016 .

Observations of neutron stars offer one way to constrain the ultra-dense EOS, with densities in the cores of neutron stars reaching $8-10\rho_{\rm sat}$. Macroscopic observations of neutron star masses Demorest2010 ; Antoniadis2013 ; Fonseca2016 and radii Ozel2009a ; Guver2010 ; Guillot2013 ; Guillot2014 ; Heinke2014 ; Nattila2016 ; Ozel2016a ; Bogdanov2016 ; Ozel2016 can be used to constrain realistic or parametric EOS, by taking advantage of the one-to-one mapping between the EOS and the mass-radius relation Lindblom1992 . Reference Ozel2010 performed the first inference of the EOS from X-ray observations. They used a sample of three radii measured from neutron stars undergoing thermonuclear bursts to infer the parameters of piecewise polytropic EOS. Reference Steiner2010 performed a similar inference for a sample extended to include three additional radii measured from transient low-mass X-ray binaries. In Ref. Steiner2013 , the sample was again extended to eight sources and a wider variety of parametric EOS were studied, including piecewise polytropic and piecewise linear models, as well as models allowing quark matter at high densities. For their baseline polytropic model, the authors constrained the radius of a $1.4M_{\odot}$ star to the range $10.4-12.9$ km. The inferred EOS constraints are recreated in Fig. 3.

New measurements over the past few years have further increased the sample of X-ray radii. In a recent analysis, Ref. Ozel2016 used a Bayesian inference scheme to combine fourteen radius measurements from X-ray bursts and quiescent low-mass X-ray binaries, while simultaneously requiring consistency with low-energy nucleon-nucleon scattering data as well as a maximum mass that is consistent with the observations of the most massive neutron stars Demorest2010 ; Antoniadis2013 ; Fonseca2016 . They inferred the parameters of a three-polytrope EOS model, which I recreate in Fig. 3, and constrained $R(1.5~{}M_{\odot})$ to lie in the range $9.9-11.2$ km Ozel2016 , slightly below the range inferred by Ref. Steiner2013 . With a more restrictive assumption of a mono-parametric EOS, i.e., assuming that all neutron stars have the same radius, Ref. Ozel2016a use a sample of 12 radius measurements to infer that the neutron star radius is $10.3\pm 0.5$ km.

While these radii measurements are starting to converge, there remain uncertainties. For a recent review of relevant systematic uncertainties and open questions, see Ref. Ozel2016 . The in-progress NASA Neutron Star Interior ExploreR (NICER) mission Gendreau2012 will provide radii measurements from pulse profile modeling of X-ray pulsars, with uncertainties potentially as small as 0.5 km. These anticipated data may provide more stringent constraints on the EOS. Additionally, a measurement of the moment of inertia from the double pulsar system, J0737$-$3039 Lyne2004 , is expected at 10% accuracy in the coming years and will offer independent constraints on the EOS Morrison2004 ; Lattimer2005 ; Kramer2009 ; Raithel2016a ; Kehl2016 .

The advent of gravitational wave observations offers a complementary avenue for constraining the EOS, independent of any of the systematics which may have biased existing measurements. Moreover, the tidal deformability represents a measurement of the quadrupolar mass distribution of the star, which has never before been measured and thus offers fresh insight into the neutron star interior. I will now turn to these constraints and insights for the remainder of the paper.

3 Tidal deformability as a probe of the neutron star radius

The tidal deformability measures the tidally-induced deformation of a neutron star, relative to the tidal potential of its binary companion. As such, it contains rich information about the stellar structure. Equation (3) shows that the dimensionless tidal deformability, $\Lambda$, can be written simply in terms of the stellar compactness and the tidal apsidal constant, which in turn depends on the compactness and the particular EOS. However, the tidal apsidal constant has been constrained to a relatively narrow range, $0.05\lesssim k_{2}^{\rm(tid)}\lesssim 0.15$, for realistic, hadronic EOS Hinderer2008 ; Hinderer2010 ; Postnikov2010 .

In a pivotal study, Ref. Yagi2013 showed that the narrow range of $k_{2}^{(\rm tid)}$ corresponds more generally to a quasi-universal relation between $\Lambda$ and the stellar compactness that holds for a variety of EOS. This relation is a corollary of the so-called “I-Love-Q” universal relations between the moment of inertia, tidal Love number, and quadrupole moment of neutron stars. Subsequently, Ref. Yagi2017 found that the $\Lambda-C$ relationship can be written, for an even larger sample of EOS, as

[TABLE]

where the coefficients were fit to be $a_{0}=0.360,a_{1}=-0.0355$, and $a_{2}=0.000705$, with errors of $\lesssim$6.5%.

Reference Annala2018 extended the work to a sample of 260,000 EOS, which they constructed to match chiral effective theory results at low densities and perturbative QCD results at high densities. In between these two regimes, the EOS were allowed to span the entire thermodynamically-consistent phase space, with either 3 or 4 polytropic segments. They found that, for a star with a fixed mass of a $1.4~{}M_{\odot}$,

[TABLE]

essentially extending the universal relation in compactness to apply more directly to radius and for a broader range of EOS. Combined with the initial gravitational wave constraint on $\Lambda(1.4~{}M_{\odot})<800$ from Ref. Abbott2017a , the authors constrained the radius of a $1.4~{}M_{\odot}$ star to be less than 13.6 km, marking one of the first constraints on the neutron star radius from GW170817.

Similar analyses have been performed based on state-of-the-art crustal models at low densities and perturbative QCD at high densities Most2018 or chiral effective field theory models Lim2018 ; Tews2018 , all of which calculated a large array of EOS models, subject to these various nuclear constraints. The ranges of tidal deformabilities predicted by the constrained EOS models were found to be consistent with the tidal deformability inferred from GW170817, but these studies point to the potential constraining power of future detections. In particular, Refs. Most2018 ; Lim2018 highlighted that a lower limit on $\widetilde{\Lambda}$ would provide strong constraints on the EOS.

In parallel to these early estimates of $R$ from the component tidal deformabilities of GW170817, it was discovered that the effective tidal deformability of the binary system can, surprisingly, also be used as a direct probe of the radius. The mass dependence that enters $\widetilde{\Lambda}$ (eq. 4) both explicitly through the masses of the component stars, as well as implicitly through $\Lambda_{1}$ and $\Lambda_{2}$, was found to be extremely weak, when the chirp mass is fixed. Because $\mathcal{M}_{c}$ is tightly constrained for GW170817, this feature renders $\widetilde{\Lambda}$ a direct probe of the radius.

Early hints of the potential of using $\widetilde{\Lambda}$ and $\mathcal{M}_{c}$ to constrain the radius can be found in Ref. Wade2014 . They showed that a dimensionful radius-like parameter can be constructed from $\widetilde{\Lambda}$, so that a measurement of $\mathcal{M}_{c}$ and $\widetilde{\Lambda}$ can be viewed as a measurement of a chirp mass and “chirp radius.” With this perspective, the inversion of $\mathcal{M}_{c}-R_{c}$ to EOS parameters becomes analogous to the traditional inverse stellar structure problem, which relates mass-radius measurements to the EOS Lindblom1992 .

Reference Raithel2018 showed the first evidence of the surprisingly tight correlation between $\widetilde{\Lambda}$ and the primary neutron star radius, which I recreate in Fig. 1. The authors found that, for a wide range of realistic EOS and a fixed chirp mass, the binary tidal deformability correlates strongly with the radius of the primary neutron star, regardless of the individual component masses. This trend can be seen in Fig. 1. In this figure, I fix the chirp mass to the central value for GW170817, $\mathcal{M}_{c}=1.186~{}M_{\odot}$ Abbott2019 . I also show, in blue, the one-dimensional posteriors of $\widetilde{\Lambda}$ digitized from the LIGO-Virgo re-analysis in Ref. Abbott2019 . Clearly, the tidal deformability from GW170817 points to a small radius of $\lesssim 13.5$ km.

Reference Raithel2018 also explored the origin of this relationship in the quasi-Newtonian regime. In particular, the authors applied a metric correction to the Newtonian tidal deformability of an $n=1$ polytropic EOS to define a “quasi-Newtonian” version of the effective tidal deformability of the binary system, $\widetilde{\Lambda}_{qN}$. Using a series expansion on $\widetilde{\Lambda}_{qN}$ in terms of the mass ratio, $q\equiv m_{2}/m_{1}$, they revealed a very weak dependence of $\widetilde{\Lambda}_{qN}$ on $q$, i.e.,

[TABLE]

where

[TABLE]

[TABLE]

and they introduced

[TABLE]

as an “effective compactness” Raithel2018 . The expression for $\widetilde{\Lambda}_{qN}$ is shown as the thin purple band in Fig. 1 for $q=0.7-1.0$, with $\mathcal{M}_{c}=1.186~{}M_{\odot}$.

Equation 7 reveals that the individual component masses enter only at $\mathcal{O}((1-q)^{2})$, while $\widetilde{\Lambda}$ depends strongly on the radius and chirp mass. Thus, for a well-constrained chirp mass, eq. (7) can be used to solve directly for the radius, with errors of at most $\sim 4\%$ for realistic neutron stars Raithel2018 .

One can further use the analytic form of eq. (7) to convert the measured posterior distribution in $\widetilde{\Lambda}$ to a distribution in radius, according to

[TABLE]

I show the corresponding distribution in Fig. 2, for a range of mass ratios. I find that the latest posteriors for $\widetilde{\Lambda}$ from GW170817 Abbott2018 imply a 90% highest posterior density interval of $9.8<R<13.2$ km (for $q=1$), with distinct likelihood peaks at $\sim$10.8 km and 12.3 km.

This range is consistent with, though slightly broader than that inferred from X-ray radii measurements for polytropic EOS (see $\S$2). Figure 2 also shows, as the gray, dashed line, the radius inferred from a sample of 12 X-ray radii measurements, under the more restrictive assumption that all neutron stars share a common radius with a mono-parametric EOS Ozel2016a . This distribution roughly agrees with the higher-significance, smaller-radius peak from GW170187. Finally, Fig. 2 indicates that the inferred radii are approximately constant across a wide range of mass ratios, as expected from the weak $q$-dependence discovered in eq. (7).

The weak-dependence on the mass ratio was also reported early on in Refs. De2018 and Zhao2018 , using a simpler approximation that employs the scaling relation, $\widetilde{\Lambda}\sim C^{-6}$. The authors of both works used the fact that $k_{2}^{(\rm tid)}$ scales roughly inversely with compactness for masses $\geq 1~{}M_{\odot}$ Hinderer2010 ; Postnikov2010 , to show that the dimensionless tidal deformability can be written as

[TABLE]

where $a\in[0.0086,0.01]$ for a large sample of polytropic EOS De2018 . The effective tidal deformability of the binary system is then

[TABLE]

where $\hat{R}$ is an effective binary radius and $f(q)$ is a weakly-dependent function of q, defined as

[TABLE]

Neglecting the $q$-dependence of $f(q)$, they then approximate

[TABLE]

with $a\prime=0.0042\pm 0.0004$. This allows a direct inversion to the radius, according to

[TABLE]

De2018 , with similar results found in Ref. Zhao2018 . I show the relationship predicted by eq. (15) in Fig. (1) as the light green band. This approximate relationship spans a broader range, but is in qualitative agreement with the more rigorous quasi-Newtonian result of Ref. Raithel2018 .

In order to obtain the most robust constraints possible from their $\widetilde{\Lambda}-R$ relationship, Ref. De2018 re-analyzed the data from GW170817, using Bayesian parameter estimation with the added requirement that both stars obey the same EOS, which was not imposed in the initial LIGO-Virgo analysis Abbott2017a . They also allowed for different priors on the mass distribution. When using the double neutron star mass distribution as the prior, they found $\widetilde{\Lambda}=245\begin{subarray}{c}+453\\ -151\end{subarray}$, corresponding, via eq. (16), to a common radius of $\hat{R}=10.9\begin{subarray}{c}+2.1\\ -1.6\end{subarray}\pm 0.2$ km. Here and in the following estimates, the 0.2 km “systematic” error comes from the uncertainty in $a\prime$ in the quasi-universal relation of eq. (15).

In an effort to make even tighter, but much more model-dependent constraints, Ref. Radice2018a performed a similar parameter estimation on GW170817 and additionally included electromagnetic data. To calculate the likelihood of the electromagnetic observations, the authors used a fitting formula from numerical relativity simulations that relates $\widetilde{\Lambda}$ to the remnant disk mass. Extant light curve models provide an estimate of the observed disk mass, which can then be compared to the value predicted by the simulated formula for a given $\widetilde{\Lambda}$. This method yielded a 90% confidence interval for $\widetilde{\Lambda}$ of (323,776) with median 487. Using eq. (16), the authors converted this to a common radius of $12.2\begin{subarray}{c}+1.0\\ -0.8\end{subarray}\pm 0.2$ km.

In a similar analysis, Ref. Coughlin2018 performed a parameter estimation using gravitational wave data, combined with radiative transfer models of the kilonova lightcurve, phenomenological fits from numerical relativity simulations of the disk and dynamical ejecta masses, and fits between the GRB jet energy and disk mass. They found $\widetilde{\Lambda}\in[292,822]$ at the 90% confidence level, with a common radius inferred via eq. (16) of $R\in[11.1,13.4]$ km, also at the 90% confidence level, with an additional 0.2 km systematic uncertainty. While these results are in rough agreement with those of Refs. Raithel2018 ; De2018 ; Zhao2018 ; Radice2018a , they rely on significant assumptions about the electromagnetic model fitting.

The $\widetilde{\Lambda}-R$ relationship was initially only tested with relatively simple, hadronic EOS. However, it remains possible that neutron-star matter undergoes a phase transition to deconfined quark matter at high densities. Depending on the exact density of the transition, a merger could could contain two neutron stars, two hybrid stars (with nuclear matter at low densities and quark matter in their cores), or one star of each. The latter case, in particular, poses a potential complication for the $\widetilde{\Lambda}-R$ relationship.

Reference Han2018 explored the effect of such a first-order phase transition on the tidal deformability, using hybrid EOS constructed with a constant-sound-speed parametrization. The authors found a much stronger dependence of $\widetilde{\Lambda}$ on $q$, in the presence of a strong phase transition. In fact, if such a phase transition indeed were present in the EOS, Ref. Han2018 showed that it may be possible to observe different values for $\widetilde{\Lambda}$ for mergers with identical chirp masses but different mass ratios. Reference Burgio2018 similarly showed that the monotonic relationship between $R(1.5~{}M_{\odot}$) and $\widetilde{\Lambda}$ is violated if there is a strong phase transition in the EOS. This is an important caveat to consider when interpreting the radii presented in this section.

4 Constraints on the maximum mass

The maximum mass of a non-rotating neutron star, $M_{\rm max}$, offers additional constraints on the EOS. Whereas the radii discussed in $\S$3 are set by the pressure near $1-2\rho_{\rm sat}$ Lattimer2001 , the maximum mass is determined by the pressure at high densities (e.g., near $\sim 8\rho_{\rm sat}$ for certain EOS Ozel2009 ) rendering it a particularly useful probe of pressures in the high-density regime.

While there exist restrictive lower limits on $M_{\rm max}$ from the observations of two nearly $2~{}M_{\odot}$ neutron stars Demorest2010 ; Antoniadis2013 ; Fonseca2016 , the upper bounds on $M_{\rm max}$ are much less constraining. A theoretical upper limit can be calculated by assuming a known EOS at low densities and forcing the sound speed to be causal at high densities, as was first demonstrated in Ref. Rhoades1974 . Setting the transition density at $1.7\rho_{\rm sat}$ results in $M_{\rm max}\sim 3.2~{}M_{\odot}$ Rhoades1974 . With more modern EOS employed up to twice the nuclear saturation density, $M_{\rm max}\sim 2.9~{}M_{\odot}$ Kalogera1996 . Most modern EOS do not exceed these limits.

There have been attempts to set observational upper limits on $M_{\rm max}$, e.g., from analyses of short GRB formation Fryer2015 ; Lawrence2015 or from the mass distribution of binary neutron stars Antoniadis2016 ; Alsing2018 . However, these methods are relatively indirect. Observations of GW170817 and its electromagnetic counterpart offer new, independent constraints. These new upper bounds point towards $M_{\rm max}\lesssim 2.3~{}M_{\odot}$, if the remnant collapses to a black hole early on. In the case of a long-lived neutron star remnant (lifetime $\gtrsim$ days), the maximum mass may be much larger.

Following a merger, the remnant mass, $M_{\rm rem}$, will determine the object’s ultimate fate. There are three possible outcomes. If $M_{\rm rem}$ exceeds the maximum mass of a differentially rotating star, the remnant will undergo prompt collapse to a black hole. If $M_{\rm rem}>M_{\rm max}$, but is below the rotating maximum mass, the remnant will survive for a short time as a hypermassive neutron star (HMNS) supported by differential rotation, or as a supramassive neutron star (SMNS) supported by rigid-body rotation. The massive neutron star will spin down until it can no longer be supported through rotational energy, and the object will collapse to a black hole. Finally, if $M_{\rm rem}<M_{\rm max}$, the remnant neutron star may survive indefinitely.

Out of these options, Ref. Margalit2017 argued that the remnant must have been a short-lived HMNS, which collapsed to a black hole after $\gtrsim 10$ ms. A longer-lived remnant would have transferred significantly more kinetic energy to its surroundings as it spun down, which is inconsistent with the observed energies of the kilonova and GRB afterglow. By combining the total binary mass, inferred from the gravitational wave signal of GW170817, with these energetic constraints, the authors inferred $M_{\rm max}\lesssim 2.17~{}M_{\odot}$, at the 90% confidence level.

In contrast, Ref. Rezzolla2018 argued for a longer-lived remnant, that is born as a HMNS but collapses to a black hole as a SMNS after redistributing its angular momentum. The authors argued that, in order to produce the observed blue kilonova, the ejecta must have had a large electron fraction, which could only be produced by the hot, polar region of a relatively long-lived remnant. However, the quick production of the short gamma-ray burst ($\sim 1$ s post-merger), implies that the remnant did not survive indefinitely, but rather collapsed to a black hole. From these arguments, the authors inferred that the remnant collapsed as a SMNS, with a mass near the mass-shedding limit (i.e., the maximum mass that can be supported by uniform rotation). Using a quasi-universal relation between the $M_{\rm max}$ and the mass-shedding limit, they concluded $M_{\rm max}\lesssim 2.16\begin{subarray}{c}+0.17\\ -0.15\end{subarray}M_{\odot}$.

Reference Shibata2017 used numerical relativity simulations to similarly infer that the remnant must have been a longer-lived massive neutron star, but surrounded by a torus. In this scenario, the remnant collapses to a black hole within milliseconds to seconds, depending on the EOS. The authors found that strong neutrino irradiation from a relatively long-lived remnant could provide one way to avoid contamination by lanthanide elements along the line of sight, and could explain the inferred low ejecta opacity. While the presence of a long-lived remnant supports a large $M_{\rm max}$, the absence of optical counterparts from the ejecta restricts the mass from being too large. The authors thus concluded that $M_{\rm max}\sim 2.15-2.25M_{\odot}$.

Finally, Ref. Ruiz2018 used GRMHD simulations to argue that the observed GRB can be explained by a HMNS that undergoes delayed collapse. In this scenario, the remnant mass can be supported by differential rotation as a HMNS, but its mass exceeds that which can be supported through uniform rotation. Using a set of quasi-universal relations inferred from GRMHD simulations for the ratio of the threshold mass for prompt collapse of a HMNS to $M_{\rm max}$ and the ratio of the mass for delayed collapse to $M_{\rm max}$, the authors constrained $M_{\rm max}\in(2.16-2.28)~{}M_{\odot}$.

These estimates of the maximum mass all point towards a common value of $M_{\rm max}\lesssim 2.3~{}M_{\odot}$, based, in part, on the timescale of the collapse of the remnant to a black hole. However, some studies have recently argued that the late-time emission from AT2017gfo can be better explained by energy injection from a long-lived neutron star remnant (with a lifetime $\gtrsim$ days) Ai2018 ; Geng2018 ; Li2018 ; Yu2018 . This scenario implies a very stiff EOS, in order to support a maximum mass $M_{\rm max}\gtrsim 2.6~{}M_{\odot}$ Yu2018 . Thus, depending on the interpretation of the electromagnetic counterpart, the implications for $M_{\rm max}$ may vary significantly. Continued observations of the lightcurve or of the polarization evolution may help resolve these differing interpretations for GW170817 Geng2018 .

In future observations, post-merger gravitational waves from the remnant may provide a clear indication of the remnant’s fate. While no such emission was detected for GW170817, it may be detectable for similar events when advanced detectors reach design sensitivity or next-generation detectors come online Abbott2017d .

5 Bayesian inferences of EOS parameters from GW170817

I have so far focused on using the neutron star radius (via $\widetilde{\Lambda}$) to constrain the pressure near $1-2~{}\rho_{\rm sat}$ or using $M_{\rm max}$ to constrain the pressure at high densities. In this section, I now turn to the third and final approach that I will discuss: using the strain data from GW170817 to constrain the EOS of neutron stars directly, via Bayesian inference schemes.

Before going further, it is useful to note why such inference schemes are necessary. Were one to have complete knowledge of the observable functions $\Lambda(\rho_{c})$ and $m(\rho_{c})$, parameterized arbitrarily by the central density of each star, this would uniquely determine the EOS and the functions could be mathematically inverted to go from $\Lambda-m$ space to pressure-density space. This is a standard example of the inverse stellar structure problem Lindblom2014 , but only applies for single stars. Recently, Ref. Lindblom2018 showed that the single-star inverse structure problem can be generalized to apply to binary systems, and that the set of observables $(m_{1},m_{2},\widetilde{\Lambda})$ can also be inverted to recover the EOS. However, the inversion is successful only in the limit of a large number of observations that span the complete $(m_{1},m_{2},\widetilde{\Lambda})$-parameter space. Bayesian inference schemes offer a way to constrain the EOS in the sparse-data limit, in which there are fewer observations than parameters describing the EOS, by incorporating Bayesian priors and even combining independent types of observations.

In Bayesian parameter estimation for a neutron star merger, the parameters are typically decomposed into those that represent a point-particle waveform, $\theta_{pp}$, and those that depend on the EOS, $\theta_{\rm EOS}$ Abbott2018 . The latter set of parameters can be represented, for example, by the binary effective tidal deformability and its correction term, $\widetilde{\Lambda}$ and $\delta\widetilde{\Lambda}$, or by the tidal deformabilities of the component stars, $\Lambda_{1}$ and $\Lambda_{2}$. Alternatively, $\theta_{\rm EOS}$ can be represented by a parameterized EOS, such as those that use piecewise polytropes or linear segments Read2009 ; Ozel2009 ; Hebeler2010 ; Steiner2010 ; Raithel2016 , spectral expansions of the adiabatic index Lindblom2010 , or constant sound-speed parameterizations Alford2013 .333For a recent discussion of the effectiveness of piecewise polytropes compared to a spectral expansion parameterization in gravitational wave inference, see Ref. Carney2018 . For a comparison with respect to a sound-speed parametrization, see Ref. Greif2019 . Both these studies highlight the problems that discontinuities in a piecewise-polytropic EOS can pose. Whereas $\S$3 focused on inferences that sampled $\widetilde{\Lambda}$ and related the results to EOS constraints, in this section, I will review Bayesian schemes that sample the EOS function directly.

Prior to the observation of GW170817, the inference techniques for connecting observables to the EOS had been well developed. In one early study, Ref. Read2009a used numerical simulations to show that the difference between a point-particle inspiral and the inspiral of two neutron stars with realistic EOS could be detected at an effective distance of 100 Mpc with Advanced LIGO. They found that, at that distance, EOS predicting neutron stars that differ in radius by 1 km could be be distinguished, and the pressure at $1.85\rho_{\rm sat}$ could be constrained to within $5\times 10^{33}$ dyn/cm2 Read2009a . With an improved waveform model calibrated to results from numerical simulations, it becomes possible to distinguish EOS predicting radii that differ by $\sim 10\%$ at distances of up to 300 Mpc Read2013 . Further numerical simulations have found that tidally-corrected effective-one-body models can accurately measure tidal effects in the early inspiral phase as well Baiotti2010 ; Baiotti2011 ; Bernuzzi2012 ; Damour2012 ; Hotokezaka2013 .

While many of these studies used Fisher matrices to determine detectability of tidal effects on the waveform Read2009a ; Read2013 ; Flanagan2008 ; Hinderer2010 ; Damour2012 , the Fisher matrix method relies on high signal-to-noise, which may be difficult to attain for gravitational wave events, even with multiple detections of neutron star-neutron star mergers Cokelaer2008 ; Vallisneri2008 . In the first fully Bayesian parameter estimation of a synthetic gravitational waveform, Ref. Del-Pozzo2013 showed that a few tens of detections will provide strong constraints on existing EOS models, by fitting $\theta_{\rm EOS}$ for the tidal deformability. Shortly thereafter, Ref. Lackey2015 performed a Bayesian inference on synthetic data that directly sampled the parameters of a piecewise polytropic EOS. They showed that the EOS pressures could be constrained to within a factor of 2 at supranuclear densities, with just one year of observing with a three-detector network, assuming 40 events per year and signal-to-noise $>8$ in each detector. Moreover, they found that most of the constraining power would come from the loudest $\sim$5 events Lackey2015 , yielding an optimistic prediction for Advanced LIGO operations. In a more recent but similar mock analysis, Ref. Abdelsalhin2018 showed that a mere three detections of neutron star mergers, observed by a four-detector network (two LIGO, one Virgo, and one KAGRA detector), would lead to the determination of the central pressures of neutron stars to within 10%, with potential errors in the pressure at $1.85\rho_{\rm sat}$ as low as 1%, for optimal scenarios.

In the LIGO-Virgo Collaboration inference of the neutron star radius from GW170817, Ref. Abbott2018 performed the Bayesian parameter estimation in both ways: by sampling the tidal deformability parameters and, in a separate analysis, by directly sampling the parameters of the EOS function. In the former case, a set of quasi-universal relations were used to link the tidal deformabilities to neutron star radii, and hence the EOS. For the case in which the EOS function was sampled directly, the authors used a spectral parameterization Lindblom2010 , and incorporated astrophysical priors on the neutron star maximum mass to further constrain the parameter space. At the 90% confidence level, they inferred radii of the two stars to be $R_{1}=10.8\begin{subarray}{c}+2.0\\ -1.7\end{subarray}$ km and $R_{2}=10.7\begin{subarray}{c}+2.1\\ -1.5\end{subarray}$ km using the tidal deformability sampling, and $R_{1}=11.9\begin{subarray}{c}+1.4\\ -1.4\end{subarray}$ km and $R_{2}=11.9\begin{subarray}{c}+1.4\\ -1.4\end{subarray}$ km for the EOS sampling Abbott2018 . The latter radii are larger as a result of the extra constraint to support a 1.97 $M_{\odot}$ neutron star, which disfavors very soft EOS. Both sets of radii are consistent with the current range of radii of 9.9-11.2 km, inferred by Ref. Ozel2016 from the composite set of X-ray radius measurements, low energy nucleon-nucleon scattering data, and the maximum mass constraint (see $\S$2).

Using the spectral EOS sampling, the LIGO-Virgo collaboration also reported direct constraints on the pressure, finding the pressure at twice the nuclear saturation density to be $3.5\begin{subarray}{c}+2.7\\ -1.7\end{subarray}\times 10^{34}$ dyn/cm2, and at $6\rho_{\rm sat}$ $9.0\begin{subarray}{c}+7.9\\ -2.6\end{subarray}\times 10^{35}$ dyn/cm2, at the 90% confidence level Abbott2018 .

In contrast to these analyses, Ref. Landry2018 argued for a non-parametric inference of the EOS. In this scheme, non-parametric priors are constructed using Gaussian processes that have been conditioned on a set of realistic EOS. In a Monte Carlo integration, synthetic EOS are drawn from these non-parametric priors and used to predict the tidal deformabilities for given component masses. These macroscopic properties are then compared with the observed data to calculate the EOS likelihoods. The authors argued that this method will more faithfully recreate complicated EOS, which may have phase transitions or other discontinuities. For GW170817, they found the pressure at twice the nuclear saturation density to be $1.35\begin{subarray}{c}+1.8\\ -1.2\end{subarray}\times 10^{34}$ dyn/cm2, for the case in which the prior was only loosely trained on the realistic EOS, and $4.73\begin{subarray}{c}+1.4\\ -2.5\end{subarray}\times 10^{34}$ dyn/cm2, for the case in which the prior was tightly trained on the sample EOS. For GW170817, this method thus produces results consistent with the LIGO-Virgo parametric analysis Abbott2018 .

I summarize all of these results in Fig. 3. The orange and blue bands represent the EOS inferred from X-ray measurements of neutron star radii in Refs. Steiner2013 and Ozel2016 , respectively (see $\S$2 for details of those analyses). Atop these existing constraints, I show the new constraints on the pressures at $2\rho_{\rm sat}$ and $6\rho_{\rm sat}$ from the parametric Bayesian inference of Ref. Abbott2018 as well as the non-parametric inference of Ref. Landry2018 . While these analyses all report constraints at different confidence levels, it is clear that the two families of observations give approximately consistent results. With additional data– both of new gravitational wave sources and of new X-ray radii with potentially smaller uncertainties, as are expected from NICER – the constraints will become even tighter.

6 Finite-temperature effects

The EOS constraints discussed in this paper so far have all corresponded to the cold EOS, for which thermal effects are ignored. However, in dynamical phenomena such as neutron star mergers, the temperature may range from this “cold” regime, to temperatures of up to 10-100 MeV Oechslin2007 . The additional thermal support provided at these temperatures may affect post-merger gravitational wave frequencies, the lifetime of the remnant, post-collapse accretion disk masses, or other observable features (see, e.g., Ref. Baiotti2009 ).

If we wish to constrain the cold EOS from observations of these features, we need to be able to disentangle the cold and thermal pressures. A small number of finite-temperature EOS have been formulated to date, which self-consistently calculate thermal effects. Among these are the LS compressible liquid drop model with a Skyrme nuclear force Lattimer1991 , the Shen relativistic mean field (RMF) thoery model with a Thomas-Fermi approximation Shen1998a , and the HS statistical model that has since been applied to $\sim$10 new RMF models and/or mass tables Hempel2010 . For a recent review on finite-temperature EOS, see Ref. Lattimer2016 . However, compared to the diversity of cold EOS models described in $\S$2, the range of physics probed by finite-temperature EOS is relatively limited.

In order to span a wider range of possible physics, many authors have employed so-called “hybrid” EOS, in which the thermal pressure of an ideal fluid is added to a cold EOS to account for heating during the merger Janka1993 . This approach is computationally simple to implement and allows the huge diversity of parametric EOS to be used in representing the cold part of the EOS. However, by representing the thermal component as an ideal fluid, hybrid EOS neglect the significant effects of degeneracy, which become especially important at high densities (for the effects of degeneracy on the thermal pressure, see, e.g., Refs. Constantinou2015 ; Constantinou2015a ; Constantinou2017a ).

Simplifying the thermal effects in this way has important consequences for the inferred observables of mergers. For example, Ref. Bauswein2010 showed that using hybrid EOS in numerical simulations of a neutron star-neutron star merger can result in post-merger gravitational wave frequencies that differ by 50-250 Hz from those predicted by more realistic EOS. The lifetime of the hypermassive remnant and post-accretion disk mass following the collapse can also differ significantly.

If progress is to be made in understanding observables from mergers, it is imperative to develop improved techniques for modeling thermal effects. Ideally, such models will allow a wide range of underlying physics in the cold EOS, and will also go beyond the ideal-fluid approximation of hybrid EOS. A first step towards this goal has recently been undertaken in Ref. Raithel2019 , in which the authors create a thermal extension that can be added to any cold EOS. They find that degenerate thermal effects can be reasonably approximated with two additional parameters beyond the ideal-fluid model, with errors of $\lesssim 30\%$ compared to more realistic calculations. While these results are promising, the full impact of these improved models on merger observables still needs to be explored.

7 Conclusions

The discovery of GW170817 has opened up a new era for constraining the neutron star EOS. In particular, the measurement of the tidal deformability of the binary system has unlocked new constraints that map directly to stellar parameters. In this paper, I have reviewed the surprising correlation between $\widetilde{\Lambda}$ and $R$. Using the re-analysis of GW170187 by the LIGO-Virgo collaboration Abbott2019 , I have shown that the implied radius of the primary neutron star is $9.8<R<13.2$ km, for $q=1$, which is consistent with measurements of the radii from X-ray observations Steiner2013 ; Ozel2016 . While the $\widetilde{\Lambda}-R$ relationship is a promising tool for hadronic EOS, it seems to break down in the event of a first-order phase transition. GW170817 also inspired new upper limits on the maximum mass of neutron stars, with several independent studies pointing towards $M_{\rm max}\lesssim 2.3~{}M_{\odot}$ and the most restrictive constraint finding, at the 90% confidence level, $M_{\rm max}\lesssim 2.17~{}M_{\odot}$ Margalit2017 , under the assumption that the remnant collapses to a neutron star within a few seconds post-merger. These analyses provide more direct observational constraints on the maximum $M_{\rm max}$ than have previously been possible, and can be used to place corresponding upper limits on the EOS pressure at high densities. Additionally, Bayesian parameter estimates allow for direct constraints on EOS, using both parametric and non-parametric methods. Finally, I reviewed the importance of improving finite-temperature EOS models, in order to be able to disentangle thermal effects from the cold EOS in interpreting future observations of neutron star mergers.

Acknowledgements.

I gratefully acknowledge useful conversations with and feedback from Feryal Özel and Dimitrios Psaltis. This work is supported by NSF Graduate Research Fellowship Program Grant DGE-1143953 and support from NASA grant NNX16AC56G.