Constraining nuclear matter parameters with GW170817

TL;DR

This paper refines constraints on nuclear matter parameters using gravitational wave data from GW170817, accounting for broader models, correlations, and uncertainties, leading to more conservative bounds on key nuclear parameters.

Contribution

It extends previous analyses by incorporating a wider range of equations of state, correlations with measured tidal deformability, and uncertainties, providing more robust bounds on nuclear matter parameters.

Findings

GW170817 bounds on $K_0$, $M_0$, and $K_{sym,0}$ are more conservative.

The analysis accounts for equation of state variations and binary mass ratio effects.

Updated tidal deformability measurement leads to refined parameter constraints.

Abstract

The tidal measurement of gravitational waves from the binary neutron star merger event GW170817 allows us to probe nuclear physics that suffers less from astrophysical systematics compared to neutron star radius measurements with electromagnetic wave observations. A recent work found strong correlation among neutron-star tidal deformabilities and certain combinations of nuclear parameters associated with the equation of state. These relations were then used to derive bounds on such parameters from GW170817 assuming that the relations and neutron star masses are known exactly. Here, we expand on this important work by taking into account a few new considerations: (1) a broader class of equations of state; (2) correlations with the mass-weighted tidal deformability that was directly measured with GW170817; (3) how the relations depend on the binary mass ratio; (4) the uncertainty from equation of state variation in the correlation relations; (5) adopting the updated tidal deformability measurement from GW170817. Upon these new considerations, we find GW170817 bounds on nuclear parameters (the incompressibility $K_0$, its slope $M_0$ and the curvature of symmetry energy $K_{\mathrm{sym},0}$ at nuclear saturation density) to be 81 MeV $\leq K_0 \leq$ 362 MeV, 1556 MeV $\leq M_0 \leq$ 4971 MeV, and -254 MeV $\leq K_\mathrm{sym,0} \leq$ 27 MeV, which are more conservative than previously found with systematic errors more properly taken into account.

Tables

Table 1. Table 1: GW170817 constraints on the incompressibility K 0 subscript 𝐾 0 K_{0} (top row), its slope M 0 subscript 𝑀 0 M_{0} (middle row), and the symmetry energy curvature K sym , 0 subscript 𝐾 sym 0 K_{\text{sym},0} for 4 different sets of priors on L 0 subscript 𝐿 0 L_{0}   Lattimer and Lim ( 2013 ); Oertel  et al. ( 2017 ); Lattimer and Steiner ( 2014 ); Tews  et al. ( 2017 ) , and Λ ~ ~ Λ \tilde{\Lambda}   Coughlin  et al. ( 2019 ); Abbott  et al. ( 2019a ) . These quantities are computed using two different methods: (i) a simple linear regression estimation described in Sec.  V.1 (labeled “Method 1” on the right column), and (ii) a comprehensive computation of the nuclear parameter posterior probability distributions described in Sec.  V.2 (labeled “Method 2” on the left column). The two methods show moderate agreement, although the first method can be seen to over-estimate the errors – thus we recommend the use of the more accurate distributions computed in method 2, which properly take into account the covariances between the parameters, as well as utilizes the full posterior distribution on Λ ~ ~ Λ \tilde{\Lambda} derived by the LIGO Collaboration  Abbott  et al. ( 2019b ) . The bounds on nuclear parameter M 0 subscript 𝑀 0 M_{0} and K sym , 0 subscript 𝐾 sym 0 K_{\text{sym},0} are weaker but more reliable than those found in  Malik  et al. ( 2018 ) due to the inclusion of scatter uncertainty in our linear regressions. The bounds on K 0 subscript 𝐾 0 K_{0} and M 0 subscript 𝑀 0 M_{0} should be taken as a rough estimate as the correlation in universal relations that were used to derive them are not large, and thus, may be subject to change with inclusion of further EoSs.

	Method 2	Method 1
	$\tilde{\Lambda }$ Posterior Distribution Abbott et al. (2019b)	70–720 Abbott et al. (2019a)	279–822 Coughlin et al. (2019)
40–62 Lattimer and Lim (2013); Lattimer and Steiner (2014); Tews et al. (2017)	69 MeV $\le K_0\le$ 352 MeV 1371 MeV $\le M_0\le$ 4808 MeV -285 MeV $\le K_{\mathrm{sym},0}\le$ 7 MeV	100 MeV $\le K_0\le$ 375 MeV 1538 MeV $\le M_0\le$ 5433 MeV -358 MeV $\le K_{\mathrm{sym},0}\le$ 23 MeV	118 MeV $\le K_0\le$ 388 MeV 1849 MeV $\le M_0\le$ 5609 MeV -298 MeV $\le K_{\mathrm{sym},0}\le$ 54 MeV
30–86 Oertel et al. (2017)	123 MeV $\le K_0\le$ 330 MeV 1884 MeV $\le M_0\le$ 4635 MeV -285 MeV $\le K_{\mathrm{sym},0}\le$ 7 MeV	45 MeV $\le K_0\le$ 398 MeV 955 MeV $\le M_0\le$ 5675 MeV -358 MeV $\le K_{\mathrm{sym},0}\le$ 23 MeV	63 MeV $\le K_0\le$ 411 MeV 1266 MeV $\le M_0\le$ 5852 MeV -298 MeV $\le K_{\mathrm{sym},0}\le$ 54 MeV

Table 2. Table 2: Similar to Table  1 , when excluding PEs, and only considering the “first method” of computing nuclear parameter constraints. Observe how the bounds upon removal of PEs show drastic improvement - showing closer agreement with Ref.  Malik  et al. ( 2018 ) (in addition to uncertainty from EoS variation), and also highlighting the effect of utilizing a large set of additional EoSs. The exception is high order nuclear parameter K sym , 0 subscript 𝐾 sym 0 K_{\text{sym},0} - showing weakened constraints due to the inclusion of uncertainty in L 0 subscript 𝐿 0 L_{0} .

	70-720 Abbott et al. (2019a)	279–822 Coughlin et al. (2019)
40–62 Lattimer and Lim (2013); Lattimer and Steiner (2014); Tews et al. (2017)	161 MeV $\le K_0\le$ 309 MeV 1506 MeV $\le M_0\le$ 3506 MeV -327 MeV $\le K_{\mathrm{sym},0}\le$ 140 MeV	182 MeV $\le K_0\le$ 324 MeV 1851 MeV $\le M_0\le$ 3723 MeV -246 MeV $\le K_{\mathrm{sym},0}\le$ 190 MeV
30–86 Oertel et al. (2017)	134 MeV $\le K_0\le$ 320 MeV 1131 MeV $\le M_0\le$ 3662 MeV -394 MeV $\le K_{\mathrm{sym},0}\le$ 168 MeV	155 MeV $\le K_0\le$ 335 MeV 1476 MeV $\le M_0\le$ 3880 MeV -313 MeV $\le K_{\mathrm{sym},0}\le$ 218 MeV

Equations

$$Q_{ij}=-\lambda {\mathcal{E}}_{ij}.$$

$$\Lambda \equiv \frac{\lambda }{M^5},$$

$$\begin{array}{rl}\Lambda & =\frac{16}{15}(1-2\bar{C}{)}^2[2+2\bar{C}(y_{\text{mbox}R}-1)-y_{\text{mbox}R}]\\ & \times \{2\bar{C}[6-3y_{\text{mbox}R}+3\bar{C}(5y_{\text{mbox}R}-8)]\\ & +4{\bar{C}}^3[13-11y_{\text{mbox}R}+\bar{C}(3y_{\text{mbox}R}-2)+2{\bar{C}}^2(1+y_{\text{mbox}R})]\\ & +3(1-2\bar{C}{)}^2[2-y_{\text{mbox}R}+2\bar{C}(y_{\text{mbox}R}-1)]\ln (1-2\bar{C}){\}}^{-1}.\end{array}$$

$$\displaystyle\begin{split}&h^{\prime\prime}+\Big{\{}\frac{2}{r}+\Big{[}\frac{2m}{r^{2}}+4\pi r(p-\epsilon)\Big{]}e^{\lambda}\Big{\}}h^{\prime}\\ &+\Big{\{}4\pi\Big{[}5\epsilon+9p+(p+\epsilon)\frac{d\epsilon}{dp}\Big{]}e^{\lambda}-\frac{6}{r^{2}}e^{\lambda}-\Big{(}\frac{d\nu}{dr}\Big{)}^{2}\Big{\}}h=0,\end{split}$$

$$\tilde{\Lambda }=\frac{16}{13}\frac{(1+12q){\Lambda }_1+(12+q)q^4{\Lambda }_2}{(1+q{)}^5},$$

$$e(n,\delta )=e(n,0)+S_2(n){\delta }^2+\mathcal{O}({\delta }^4),$$

$$\begin{array}{rl}e(n,0)& =e_0+\frac{K_0}{2}{y}^2+\frac{Q_0}{6}{y}^3+\mathcal{O}(y^4),\\ S_2(n)& =J_0+L_{0}y+\frac{K_{\mathrm{sym},0}}{2}{y}^2+\mathcal{O}(y^3),\end{array}$$

$$M_0$$

$$C(x,y)=\frac{{\sigma }_{xy}}{\sqrt{{\sigma }_{xx}{\sigma }_{yy}}},$$

$$\sigma_{xy}=\frac{1}{N}\sum\limits^{N}_{i=0}x_{i}y_{i}-\frac{1}{N^{2}}\Big{(}\sum\limits^{N}_{i=0}x_{i}\Big{)}\Big{(}\sum\limits^{N}_{i=0}y_{i}\Big{)}.$$

$$\frac{K_0}{\text{MeV}}+\alpha \frac{L_0}{\text{MeV}}=0.1086_{-0.02064}^{+0.02172}\tilde{\Lambda }+299.1_{-64.60}^{+72.97},$$

$$\frac{M_0}{\text{MeV}}+\beta \frac{L_0}{\text{MeV}}=1.488_{-0.2038}^{+0.2456}\tilde{\Lambda }+3929_{-990.2}^{+1226},$$

$$\frac{K_{\text{sym},0}}{\text{MeV}}+\gamma \frac{L_0}{\text{MeV}}=0.2915_{-0.004080}^{+0.007287}\tilde{\Lambda }-259.1_{-118.9}^{+67.36}.$$

$$P(\tilde{\Lambda },K_{\text{sym},0})=\frac{1}{2\pi \sqrt{|\Sigma |}}{e}^{-\frac{1}{2}(x-\mu {)}^T{\Sigma }^{-1}(x-\mu )},$$

$$P(Y|{\tilde{\Lambda }}_{\text{obs}})\sim \mathcal{N}\left({\mu }_Y+\frac{{\sigma }_Y}{{\sigma }_{{\tilde{\Lambda }}_{\text{obs}}}}C({\tilde{\Lambda }}_{\text{obs}}-{\mu }_{{\tilde{\Lambda }}_{\text{obs}}}),(1-C^2){\sigma }_Y^2\right).$$

$$P(K_0+\alpha L_0)=\underset{-\infty }{\overset{\infty }{\int }}P(K_0+\alpha L_0|\tilde{\Lambda })P_{\text{LIGO}}(\tilde{\Lambda })d\tilde{\Lambda },$$

$$P(K_0)=\underset{-\infty }{\overset{\infty }{\int }}P(K_0+\alpha L_0)P(L_0)d{L}_0,$$

$$P(\epsilon)=\left\{\begin{array}[]{ll}P_{\text{tr}}&(\epsilon_{1}\leq\epsilon\leq\epsilon_{2})\\ P_{\text{tr}}+c_{\text{s}}^{2}(\epsilon-\epsilon_{2})&(\epsilon>\epsilon_{2})\end{array}\right.$$

Full text

I Introduction

One of the largest mysteries in nuclear physics comes from the determination of the equation of state (EoS) of ultra-dense nuclear matter, found exclusively in neutron stars (NSs). Many useful relations, such as the one between mass and radius, depend strongly on the EoS, and are vital to the study of nuclear physics to constrain EoSs for supranuclear matter and model-independent parameters that characterize such EoSs. Indeed, the mass-radius measurement of NSs via X-ray observations have been used to obtain constraints on nuclear matter EoSs Guver and Ozel (2013); Ozel et al. (2010); Steiner et al. (2010); Lattimer and Steiner (2014); Ozel and Freire (2016).

Recently, gravitational waves (GWs) from a binary NS merger have been detected (GW170817) Abbott et al. (2017a), which can also be used to probe nuclear physics Abbott et al. (2019a, 2018); Paschalidis et al. (2018); Burgio et al. (2018); Malik et al. (2018). This is mainly because as two NSs in a binary system inspiral due to GW emission, each of them become tidally deformed in response to the tidal gravitational field created by the companion. Such a tidal effect is characterized by the tidal deformability Flanagan and Hinderer (2008) which depends strongly on the underlying EoSs. In fact, the leading tidal parameter entering in the gravitational waveform is given by a mass-weighted combination of the two tidal deformabilities $\tilde{\Lambda}$ associated with each NS. The LIGO Scientific Collaboration and the Virgo Collaboration (LVC) recently placed a 90% credible bound on $\tilde{\Lambda}$ as $70\leq\tilde{\Lambda}\leq 720$ Abbott et al. (2019a) (see De et al. (2018) for similar bounds). Coughlin et al. Coughlin et al. (2019) further combined numerical relativity simulations with electromagnetic counterpart signals for GW170817 and derived $279\leq\tilde{\Lambda}\leq 822$ (a similar bound was also derived in Radice et al. Radice et al. (2018)). Such bounds on $\tilde{\Lambda}$ have also been mapped to those on the NS radius Annala et al. (2018); Lim and Holt (2018); Bauswein et al. (2017); De et al. (2018); Most et al. (2018).

Given that all of the EoSs proposed so far use certain approximations, one informative approach is to directly measure nuclear physics parameters which parameterize EoSs in a model-independent way. One way to obtain such a parameterization is to Taylor expand the energy per nucleon of asymmetric nuclear matter about the saturation density111Other ways of parameterizing EoSs include piecewise polytropes Read et al. (2009); Lackey and Wade (2015); Carney et al. (2018) and spectral EoSs Lindblom (2010); Lindblom and Indik (2012, 2014); Lindblom (2018); Abbott et al. (2018). See also Landry and Essick (2019) for a non-parametric inference of EoSs with GW170817.. Taylor-expanded coefficients include the symmetry energy’s slope $L_{0}$, the incompressibility $K_{0}$, its slope $M_{0}$ and the curvature of symmetry energy $K_{\mathrm{sym},0}$.

Interestingly, approximate universal relations exist among nuclear physics parameters mentioned above and NS radius at a given mass Alam et al. (2016) (see e.g. Sotani et al. (2014); Silva et al. (2016) for other universal relations involving nuclear parameters). The authors found that while individual nuclear parameters are only weakly correlated with the stellar radius, linear combinations of the form $K_{0}+\alpha L_{0}$ and $M_{0}+\beta L_{0}$ become highly correlated, where $\alpha$ and $\beta$ are chosen such that the correlation becomes maximum.

Such work was recently extended by Malik et al. Malik et al. (2018) by considering correlations with individual NS tidal deformabilities. By taking these relations to be exact and assuming individual NS masses from GW170817 to be $m_{1}=1.40M_{\odot}$ and $m_{2}=1.33M_{\odot}$, Ref. Malik et al. (2018) utilized existing measurements on tidal deformability from GW170817 Abbott et al. (2017b); Radice et al. (2018) and $L_{0}$ Abbott et al. (2019a); Oertel et al. (2017); Lattimer and Steiner (2014) to derive constraints on the nuclear incompressibility and the symmetry energies’ curvature at saturation density to be $2254\text{ MeV}\leq M_{0}\leq 3631\text{ MeV}$ and $-112\text{ MeV}\leq K_{\text{sym},0}\leq-52\text{ MeV}$, respectively.

This important first-step work of Ref. Malik et al. (2018) needs to be improved in various ways. In this paper, we propose an extension upon this work by taking into account at least the following five points of interest. First, we consider a broader class of EoSs by phenomenologically varying nuclear parameters. Second, we consider correlations among the mass-weighted tidal deformability (instead of the individual tidal deformabilities) and nuclear parameters for various mass ratios. This allows us to eliminate the need to choose specific NS masses $m_{1}$ and $m_{2}$, as was done in Ref. Malik et al. (2018). Third, instead of assuming perfect linear regression between nuclear parameters and $\tilde{\Lambda}$, the uncertainty from scatter (corresponding to the EoS variation in the approximate universal relations) is taken into account, including the covariances among parameters. Fourth, we use the recent updated posterior distribution of the dominant tidal deformability $\tilde{\Lambda}$ by LVC Abbott et al. (2019a). Finally, we investigate constraints on the incompressibility $K_{0}$ in addition to its slope $M_{0}$ and the curvature of symmetry energy $K_{\text{sym},0}$.

I.1 Executive Summary

Let us summarize important results for busy readers. First, we find new universal relations between $\tilde{\Lambda}$ and $K_{0}$, $M_{0}$, or $K_{\text{sym},0}$ (bottom panel of Fig. 1) for a number of mass ratios allowed by GW170817. Contrary to previous work, we find low-order nuclear parameters $K_{0}$ and $M_{0}$ to have very poor correlations, due to the inclusion of a broad new class of EoSs.

Additionally, we studied similar universal relations between $\tilde{\Lambda}$ and linear combinations of nuclear parameters (top panel of Fig. 1), such as $K_{0}+\alpha L_{0}$, $M_{0}+\beta L_{0}$, and $K_{\text{sym},0}+\gamma L_{0}$. We found that such relations typically have a stronger correlation than that in the case of individual nuclear parameters. This is consistent with the findings of Ref. Malik et al. (2018) on correlations between nuclear parameters and individual tidal deformabilities, though the correlations presented here are much lower than that reported in the previous work. Contrary to Ref. Malik et al. (2018) where coefficients are chosen such that correlation is maximal, we choose coefficients $\alpha=2.27$, $\beta=24.28$, and $\gamma=0$. To avoid the propagation of uncertainties from $L_{0}$, we manually choose $\alpha$ and $\beta$ to be as small as possible, while keeping in mind that the correlation with $\tilde{\Lambda}$ must be large enough to determine bounds on nuclear parameters. We arbitrarily choose $\alpha$ and $\beta$ such that correlations are $0.50$ to give one example of the derived bounds. The parameter $\gamma$ was chosen to be [math] in order to neglect the additional uncertainty accrued by the addition of $L_{0}$, possible in this case only due to the high correlations between $K_{\text{sym},0}$ and $\tilde{\Lambda}$.

Figure 2 presents 90% confidence interval on $K_{\text{sym},0}$ after GW170817, based on the universal relation in Fig. 1. In the computation of these above bounds, the posterior probability distribution on $\tilde{\Lambda}$ as derived by the LIGO Collaboration Abbott et al. (2019b) was used. In particular, we find such bounds to be -285 MeV $\leq K_{\mathrm{sym,0}}\leq$ 7 MeV 222The constraint on $K_{\text{sym,0}}$ bears a close resemblance to that in Refs. Margueron and Gulminelli (2018); Mondal et al. (2018). at the $90\%$ confidence interval. Additionally, we find bounds on $K_{0}$ and $M_{0}$ to be 69 MeV $\leq K_{0}\leq$ 352 MeV and 1371 MeV $\leq M_{0}\leq$ 4808 MeV. Such results are much weaker than the results found in Ref. Malik et al. (2018), born from the inclusion of systematic errors from a broader class of EoSs and the scatter uncertainty from EoS variation on universal relations. These results lead us to conclude that it is important to account for the large systematic errors accrued from a wider range of valid EoSs and EoS variation in the approximate universal relations.

The organization of this paper is as follows. We begin with complementary background material on NS tidal deformability in Sec. II. We continue on discussing the standard asymmetric nuclear matter parameters, and their resulting EoSs and mass-radius relations in Sec. III. We next examine the correlations between nuclear matter parameters and mass-weighted tidal deformability in Sec. IV and further use these results to derive constraints on such nuclear parameters in Sec. V. We conclude in Sec. VI by discussing our results and give possible avenues for future work. Throughout this paper, we have adopted geometric units of $G=c=1$, unless otherwise stated.

II Neutron star tidal deformability

We begin by reviewing how one can extract internal structure information of NSs via GW measurement. In the presence of a neighboring tidal field $\mathcal{E}_{ij}$, such as the binary NS system found in GW170817, NSs tidally deform away from sphericity and acquire a non-vanishing quadrupole moment $Q_{ij}$ that is characterized by the tidal deformability $\lambda$ Flanagan and Hinderer (2008); Hinderer (2008); Yagi and Yunes (2013):

[TABLE]

Such tidal deformability can be made dimensionless as:

[TABLE]

with $M$ representing the stellar mass. $\Lambda$ can be calculated via the following expression Hinderer (2008); Damour and Nagar (2009); Yagi and Yunes (2013):

[TABLE]

Here $\bar{C}\equiv M/R$ is the stellar compactness with $R$ representing the NS radius, and $y_{\mbox{\tiny R}}\equiv y(R)$ with $y(r)\equiv rh^{\prime}(r)/h(r)$, where a prime stands for taking a derivative with respect to the radial coordinate $r$. $h$ represents the quadrupolar part of the $(t,t)$ component of the metric perturbation satisfying the following differential equation:

[TABLE]

with background metric coefficients $e^{\nu}=g_{tt}$ and $e^{\lambda}=(1-2m/r)^{-1}=g_{rr}$, while $p$ and $\epsilon$ represent pressure and energy density respectively.

The above differential equation can be solved as follows. First, one needs to prepare unperturbed background solutions by choosing a specific EoS, or $p(\epsilon)$, and solve a set of Tolman-Oppenheimer-Volkoff (TOV) equations with a chosen central density (or pressure) and appropriate boundary conditions (the exterior metric being the Schwarzschild one). The stellar radius is determined from $p(R)=0$ while the mass is given by $M=m(R)=4\pi\int_{0}^{R}\epsilon(r)\,r^{2}dr$. Having such solutions at hand, one then plugs them into Eq. (4) and solves it with the boundary condition $y(0)=2$ Hinderer (2008).

Because there are two NSs in a binary, two tidal deformabilities $\Lambda_{1}$ and $\Lambda_{2}$ associated with each star enter in the gravitational waveform. However, extracting such parameters independently is challenging due to the strong correlation between them333One way to cure this problem is to use universal relations between them Yagi and Yunes (2016, 2017); De et al. (2018); Zhao and Lattimer (2018).. Thus, one can instead measure the dominant tidal parameter in the waveform, corresponding to the mass-weighted average tidal deformability given by Flanagan and Hinderer (2008):

[TABLE]

where $q\equiv m_{2}/m_{1}(<1)$ is the mass ratio between two stars.

III Nuclear matter parameters and equations of state

III.1 Asymmetric Nuclear Matter Parameters

Here we review a generic method of parameterizing EoSs. Our starting point is expanding the energy per nucleon $e$ of asymmetric nuclear matter with isospin symmetry parameter $\delta\equiv(n_{n}-n_{p})/n$ (with $n_{p}$ and $n_{n}$ representing the proton and neutron number densities respectively and $n\equiv n_{p}+n_{n}$) about $\delta=0$ (symmetric nuclear matter case) as Vidaña et al. (2009):

[TABLE]

where $e(n,0)$ corresponds to the energy of symmetric nuclear matter. $e(n,0)$ and $S_{2}(n)$ can then be characterized by once again expanding about the saturation density $n_{0}$ as:

[TABLE]

where $y\equiv(n-n_{0})/3n_{0}$. Here, the coefficients are known as the energy per particle $e_{0}$, incompressibility coefficient $K_{0}$, third derivative of symmetric matter $Q_{0}$, symmetry energy $J_{0}$, its slope $L_{0}$, and its curvature $K_{\mathrm{sym,0}}$ at saturation density, respectively. Following Refs. Alam et al. (2014); Malik et al. (2018), we further introduce the slope of the incompressibility:

[TABLE]

In this paper, we investigate correlations between the various nuclear parameters $L_{0}$, $K_{0}$ $M_{0}$, $K_{\text{sym},0}$ and the mass-weighted average tidal deformability $\tilde{\Lambda}$ in order to derive bounds on nuclear parameters from GW170817. Bounds on $M_{0}$ and $K_{\mathrm{sym,0}}$ have previously been derived in Ref. Malik et al. (2018) using GW170817, which we revisit in this paper. Current experiments and astrophysical observations place bounds on $L_{0}$ as $40\text{ MeV}<L_{0}<62\text{ MeV}$ Lattimer and Lim (2013); Lattimer and Steiner (2014); Tews et al. (2017), and $30\text{ MeV}<L_{0}<86\text{ MeV}$ Oertel et al. (2017).

III.2 Equations of State

The structure of a NS and its tidal interactions in a binary system rely heavily on the underlying EoS of nuclear matter. Because of this, we employ a wide range of 120 different nuclear models in our analysis. These EoSs can be classified into three broad categories: 24 non-relativistic EoSs with Skyrme-type interaction, 9 RMF EoSs, and 88 EoSs derived through phenomenological variation. Following Ref. Read et al. (2009), the high-density core EoSs listed above are all matched to the low-density EoS of Douchin and Haensel Douchin and Haensel (2001) at the transition density $\epsilon_{\text{tr}}$ such that the pressures are equivalent.

The EoSs in the first two classes are used also in Alam et al. (2016); Malik et al. (2018). The Skyrme models used here are: SKa, SKb Köhler (1976), SkI2, Sk13, SkI4, SkI5 Reinhard and Flocard (1995), SkI6 Nazarewicz et al. (1996), Sly230a Chabanat et al. (1997), Sly2, Sly9 Chabanat (1995), Sly4 Chabanat et al. (1998), SkMP Bennour et al. (1989), SkOp Reinhard (1999), KDE0V1 Agrawal et al. (2005), SK255, SK272 Agrawal et al. (2003), Rs Friedrich and Reinhard (1986), BSK20, BSK21 Goriely et al. (2010), BSK22, BSK23, BSK24, BSK25, BSK26 Goriely et al. (2013). On the other hand, the RMF models selected are BSR2, BSR6 Dhiman et al. (2007); Agrawal (2010), GM1 Glendenning (1991), NL3 Lalazissis et al. (1997), NL3$\omega\rho$ Carriere et al. (2003), TM1 Sugahara and Toki (1994), DD2 Typel et al. (2010), DDH$\delta$ Gaitanos et al. (2004), DDME2 Typel and Wolter (1999).

One of the new EoS classes that we consider is the phenomenological EoSs (PEs). To construct these EoSs, we followed the formalism of Ref. Margueron et al. (2018) by randomly sampling nuclear parameters $J_{0}$, $K_{0}$, $L_{0}$, $Q_{0}$ and $K_{\mathrm{sym,0}}$ as found in Table I of the above reference. Following this, nonphysical EoSs with acausal structure ($v_{s}>c$), or having decreasing pressure as a function of density were removed.

Figure 3 presents the relations among the NS mass, radius and tidal deformability for selected EoSs in different classes mentioned above. Observe that RMF EoSs tend to produce NSs with larger radii and maximum mass than those for Skyrme-types, while the PE ones generate NSs with a wide range of properties.

IV Correlations between tidal deformability and Nuclear Parameters

In this section, we study correlations among nuclear parameters and tidal deformability, where the latter can be measured from GW observations. The amount of correlation between two variables $x$ and $y$ with $N$ data points can be quantified by the Pearson correlation coefficient $C$ defined by:

[TABLE]

where the covariances $\sigma_{xy}$ are given by:

[TABLE]

$C=1$ represents absolute correlation, while $C=0$ corresponds to having no correlation.

IV.1 $\tilde{\Lambda}$ versus Nuclear Parameters

Reference Malik et al. (2018) first studied the universal relations between nuclear parameters and the tidal deformability for isolated neutron stars. The authors then map this to the GW measurement on $\tilde{\Lambda}$ by using yet another universal relation between $\tilde{\Lambda}$ and $\Lambda_{1.4}$ for a specific choice of masses in a binary neutron star that is consistent with GW170817. However, the mass ratio $q\in[0.73,1.00]$ Abbott et al. (2017a) for this event has not been measured very precisely (the lower bound of this constraint has recently been improved to $0.8$ in Ref. Coughlin et al. (2019)), and the question arises as to whether such relation holds for various $q$. As we show in Appendix A, indeed the universal relation is highly insensitive to the choice of $q$. This suggests that there are universal relations between nuclear parameters and $\tilde{\Lambda}$ for a given chirp mass $\mathcal{M}$ which has been measured with high accuracy for GW170817. Finding these universal relations is the focus of this section. Universal relations involving $\tilde{\Lambda}$ are, in some sense, practically more useful than those with $\Lambda_{1.4}$, because the former is a quantity which can be directly measured from GW observations.

Figure 4 shows the correlations between nuclear parameters ($L_{0}$, $K_{0}$, $M_{0}$, $K_{\text{sym},0}$), and the mass-weighted average tidal deformability $\tilde{\Lambda}_{q}$ evaluated at mass ratios of $q=0.73$, $0.87$ and $1.00$. The linear regression shown in each panel represents the best fit line describing the relation between nuclear parameters and $\tilde{\Lambda}$. Observe that $K_{0}$ and $M_{0}$ show very poor correlations, resulting from a disconnect between PEs and EoSs found in Ref. Malik et al. (2018). On the other hand, higher order parameter $K_{\text{sym},0}$ sees a fairly strong correlation of $\sim 0.80$. It is noted that PEs typically have values of $K_{0}$ that are much lower than those for Skyrme or RMF EoSs, while $M_{0}$ is much higher, and $L_{0}$ and $K_{\text{sym},0}$ are very similar. Let us emphasize that we have restricted to physically valid PEs which have increasing pressure, and this is why we do not have PEs with e.g. $M_{0}<2500$ MeV 444This does not mean that Skyrme and RMF EoSs with $M_{0}<2500$ are nonphysical.. The above finding indicates a necessity in using a large number of EoSs as nuclear parameters can take on a much wider range of values than considered in Malik et al. (2018). Observe also that the behavior of the scattering and the amount of correlation found in Fig. 4 is not very sensitive to $q$. This can also be seen from Fig. 5, where correlations between various nuclear parameters and $\tilde{\Lambda}$ are plotted as a function of mass ratio q.

IV.2 $\tilde{\Lambda}$ versus linear combinations of nuclear parameters

References Alam et al. (2016); Malik et al. (2018) report that correlations among nuclear parameters and NS observables become stronger if one considers certain combinations of the former, which we study here. In Refs. Oertel et al. (2017); Lattimer and Steiner (2014); Tews et al. (2017), tight constraints on the slope of the symmetry energy $L_{0}$ were derived. Thus we focus on constraining the incompressibility $K_{0}$, its slope $M_{0}$, and the symmetry energies’ curvature $K_{\text{sym},0}$, utilizing prior bounds on $L_{0}$ and $\tilde{\Lambda}$ by considering linear combinations of the form $K_{0}+\alpha L_{0}$, $M_{0}+\beta L_{0}$, and $K_{\text{sym},0}+\gamma L_{0}$ with some coefficients $\alpha$, $\beta$ and $\gamma$. In previous literature Alam et al. (2016); Malik et al. (2018), these coefficients are chosen such that correlations become maximum.

Figure 5 presents the correlations between $\tilde{\Lambda}$ and linear combinations of nuclear parameters as a function of mass ratio $q$. We found that the values of $\alpha$ and $\beta$ which give maximal correlation are unnecessarily large. For practical purposes, we choose here $\alpha=2.27$ and $\beta=24.28$, such that a correlation of 50% in the universal relations is achieved. For $\gamma$, we use $\gamma=2.63$ which maximizes the correlation, as was done previously (see Sec. V for more details). For reference, we also show correlations involving single nuclear parameters. Observe that the former correlations are much stronger than the latter (except for $K_{0}+\alpha L_{0}$ whose correlation is comparable to that of $K_{0}$) and remain to be strong over the acceptable region of mass ratio. This implies that our choice of $q$ when calculating bounds on nuclear parameters does not matter significantly. Therefore, we consider universal relations evaluated at the central mass ratio of $q=0.87$, shown in Figs. 1 and 6. Also notice how linear combinations involving high-order nuclear parameter $K_{\text{sym},0}$ continue to significantly outperform lower-order parameters.

V Constraints on nuclear matter parameters

Let us now use the approximate universal relations among combined nuclear parameters and $\tilde{\Lambda}$ to derive bounds on the former from the measurement of the latter with GW170817. In this section, we detail the process used to estimate nuclear parameter bounds, taking into account the EoS scattering uncertainty. We offer two alternative methods of accomplishing this. In Sec. V.1, we offer a crude estimation of the constraints by finding linear regressions between the nuclear parameters and $\tilde{\Lambda}$. We estimate $90\%$ confidence integrals on such regressions which allows us to predict bounds on nuclear parameters. The linear regressions provide ready-to-use type results that can easily be implemented as the measurement on $\tilde{\Lambda}$ from GW170817 are updated. In Sec. V.2, we detail a more comprehensive analysis in which we first compute the 2-dimensional probability distribution between the nuclear parameters and $\tilde{\Lambda}$. We then combine this with the probability distribution on $\tilde{\Lambda}$ computed by Ref. Abbott et al. (2019b) to estimate the posterior distribution on nuclear parameters $K_{0}$, $M_{0}$, and $K_{\text{sym},0}$.

V.1 Constraint Estimation via Linear Regressions

In this simple error analysis, we first construct linear regressions of the form $(a\pm\delta_{a}^{\pm})\tilde{\Lambda}+(b\pm\delta_{b}^{\pm})$ on the relations evaluated at the central mass ratio of $q=0.87$ with the “90%” error on the slope and $y$-intercept as follows:

[TABLE]

[TABLE]

[TABLE]

The uncertainties on the slope and $y$-intercept, $\delta_{a}^{\pm}$ and $\delta_{b}^{\pm}$, are found by varying the upper and lower error bars throughout the parameter space, selecting only combinations of $\delta_{a}^{\pm}$ and $\delta_{b}^{\pm}$ which form “90% error lines” $(a\pm\delta_{a}^{\pm})\tilde{\Lambda}+(b\pm\delta_{b}^{\pm})$ containing 90% of the data points between them. Further, we choose the “best fit” 90% error lines by minimizing the residual sum of squares, $\sum_{1=1}^{n}(y_{i}-f(x_{i}))^{2}$, as denoted by the dashed black lines in Figs. 1, and 6. For reference, the covariances $\sigma_{ab}$ from Eq. (9) between $a$ and $b$ are found to be approximately $0.7274$, $124.5$, and $0.4235$ for Eqs. (11)–(13), respectively. Using this method of uncertainty prediction, we find a 90% confidence interval on the value of $b$ and $a$, allowing us to account for the EoS scatter in the universal relations when deriving bounds on nuclear parameters from GW170817, as we will study next.

Let us now use Eqs. (11)–(13) to derive bounds on $K_{0}$, $M_{0}$, and $K_{\text{sym},0}$, as was done in Ref. Malik et al. (2018). We utilize prior bounds obtained from nuclear experiments and astrophysical observations as $L_{0}\in[40,62]$ MeV Lattimer and Lim (2013) and $L_{0}\in[30,86]$ MeV Oertel et al. (2017); Lattimer and Steiner (2014); Tews et al. (2017), as well as tidal deformability ranges of $\tilde{\Lambda}\in[70,720]$ Abbott et al. (2019a) and $\tilde{\Lambda}\in[279,822]$ Coughlin et al. (2019). Utilizing the 90% confidence interval’s range on y-intercepts, we find constraints on $K_{0}$, $M_{0}$, and $K_{\text{sym},0}$ within priors of $L_{0}$ and $\tilde{\Lambda}$ such that minimal and maximal values of nuclear parameters are obtained. Therefore, 2 constraints on $\tilde{\Lambda}$ and 2 constraints on $L_{0}$ allow us to derive 4 possible constraints on each nuclear parameter $K_{0}$, $M_{0}$, and $K_{\text{sym},0}$. This particular method of estimating the probability distribution is conservative by nature, and also takes into account the uncertainty from scatter in our relations.

The top panels of Fig. 7 show comparisons between estimated nuclear parameter limits, while the central panels show constraint ranges (maximum value minus minimum value) as the linear combination coefficient ($\alpha$, $\beta$, or $\gamma$) is increased. The bounds are stronger if the ranges are smaller. For comparison, the bottom panels display the correlation between the nuclear parameter combinations and $\tilde{\Lambda}$. Observe that the bounds become weaker as one increases the coefficients, as we are introducing an additional source of uncertainty from $L_{0}$. Does this mean that it is always better to set the coefficients to 0 and consider universal relations with individual nuclear parameters? The answer is no because correlations are too small when $\alpha=\beta=0$, as can be seen from the bottom panels of Fig. 7. If such correlations are too small, the relations can easily be affected by the addition of new EoSs and the bounds derived from these relations become unreliable.

Therefore, we need to find the balance between having large enough correlations and yet to have reasonable bounds on the nuclear parameters. Regarding $\alpha$ and $\beta$, notice that bounds on $K_{0}$ and $M_{0}$ increase approximately linearly with the coefficients, while correlations with $\tilde{\Lambda}$ quickly asymptote to values of $\sim 0.60$. Thus we choose $\alpha=2.27$ and $\beta=24.28$ such that correlations evaluated at central mass ratio $q=0.87$ are an arbitrary value of $C=0.50$, chosen to keep correlations as high as possible, while keeping $\alpha$ and $\beta$ as small as possible to avoid the propagation of uncertainty in $L_{0}$. Regarding $\gamma$, because $K_{\text{sym},0}$ starts off with strong correlation at $\gamma=0$, we choose this value to remove any additional uncertainty in $\gamma$ and $L_{0}$ from our calculations (Note this can not be done for the cases of $K_{0}$ and $M_{0}$ due to weak individual correlations with $\tilde{\Lambda}$). Observe that the coefficient choices discussed in Ref. Malik et al. (2018), to maximize correlations to the level of $0.8$ and beyond is not necessarily applicable to every situation. As seen in Fig. 7, high correlations are unobtainable for linear combinations involving $K_{0}$ and $M_{0}$, yielding no bounds under such a selection criteria. Instead, reducing the threshold to $0.50$ returns constraints as shown below, albeit being less reliable.

Table 1 summarizes the bounds on the nuclear parameters with these fiducial choices of $\alpha$, $\beta$ and $\gamma$, using both this method of constraint estimation, and the method described in Sec. V.2. The constraints on $M_{0}$ and $K_{\mathrm{sym},0}$ are additionally visualized in Fig. 2. Notice how our conservative constraints (found by using the largest-range priors on both $L_{0}$ and $\tilde{\Lambda}$) on the slope of incompressibility and the curvature, 955 MeV $\leq M_{0}\leq$ 5675 MeV and -358 MeV $\leq K_{\mathrm{sym,0}}\leq$ 23 MeV, are much weaker than those found in Ref. Malik et al. (2018) (see Fig. 2), due to the consideration of EoS scatter uncertainty, and of additional PEs with a wider range of nuclear values. We observe that the constraints derived here on $K_{\text{sym,0}}$ show good agreement with that of Refs. Margueron and Gulminelli (2018); Mondal et al. (2018). Let us emphasize that the bounds on $K_{0}$ and $M_{0}$ should be considered as rough estimates, as the correlation of 0.50 is not very large; thus these bounds are more easily affected by inclusion of yet additional EoSs than the bounds on $K_{\mathrm{sym,0}}$.

V.2 Constraint Estimation via LIGO Posterior Distributions

In this section, we offer a more comprehensive method of estimating nuclear matter constraints than was found in Sec. V.1. Previously, a rough estimate on the nuclear matter constraints was computed by finding linear regressions between $\tilde{\Lambda}$ and nuclear parameters. By estimating the $90\%$ errors on these lines, bounds on the nuclear parameters were manually approximated. In this section, we improve upon this method by (i) properly taking into account the covariance between $\tilde{\Lambda}$ and nuclear parameters by generating a multivariate probability distribution, and (ii) taking into account the full posterior probability distribution on $\tilde{\Lambda}$ as derived by the LIGO Collaboration Abbott et al. (2019b).

We begin by generating the 2-dimensional probability distribution between $\tilde{\Lambda}$ and the nuclear parameters, taking into account the specific covariances between them. For the example of $K_{\text{sym},0}$ the distribution is given by:

[TABLE]

where $x$ is the 2-dimensional vector containing $\tilde{\Lambda}$ and the given nuclear parameter, $\mu$ is the 2-dimensional vector containing the expected values of $x$, and $\Sigma$ is the $2\times 2$ covariance matrix defined with elements given by Eq. (10). This distribution is displayed in Fig. 8 for each nuclear parameter. Notice here the high degree of covariance between the variables used in this analysis - indicative of the importance for using this method of constraint extraction.

Following this, we compute the conditional probability distributions on nuclear matter parameters given a tidal observation of $\tilde{\Lambda}_{\text{obs}}$. Following Ref. Jensen (2007), the one-dimensional conditional probability distribution on nuclear parameter $Y$ is then given by

[TABLE]

In the above expression, $\mathcal{N}(\mu,\sigma^{2})$ is the normal distribution with mean and variance $\mu$ and $\sigma^{2}$, and $\mu_{\mbox{\tiny A}}$ and $\sigma_{\mbox{\tiny A}}^{2}$ are the mean and variances of $Y$ and $\tilde{\Lambda}_{\text{obs}}$.

Next, we extract the one-dimensional probability distributions on $K_{0}+\alpha L_{0}$, $M_{0}+\beta L_{0}$, and $K_{\text{sym},0}$ by combining the one-dimensional conditional distributions $P(Y|\tilde{\Lambda})$ found in Eq. (15) with the probability distribution $P_{\text{LIGO}}(\tilde{\Lambda})$ on $\tilde{\Lambda}$ derived by the LIGO Collaboration in Ref. Abbott et al. (2019b) for GW170817, shown in Fig. 9. For example, the posterior probability distribution on $K_{0}+\alpha L_{0}$ is given by:

[TABLE]

and similarly for $M_{0}+\beta L_{0}$ and $K_{\text{sym},0}$. Additionally, to find the probability distributions on $K_{0}$ and $M_{0}$, we perform one last integration over the prior probability distribution of $L_{0}$, assumed to be Gaussian with standard deviation $\sigma=\frac{1}{2}(80+36)$ and mean $\mu=\frac{1}{2}(80-36)$ Oertel et al. (2017) (or $\sigma=\frac{1}{2}(62+40)$ and $\mu=\frac{1}{2}(62-40)$ Lattimer and Lim (2013); Lattimer and Steiner (2014); Tews et al. (2017) for the alternative priors on $L_{0}$). For example, the probability distribution on $K_{0}$ is given by:

[TABLE]

with $\alpha=2.27$.

The results of these computations are shown in Fig. 10 for the more conservative priors on $L_{0}$. We observe that $K_{0}$, $M_{0}$, and $K_{\text{sym},0}$ now obey distributions that look like skewed Gaussians centered at $K_{0}=208^{+86}_{-85}\text{ MeV}$, $M_{0}=3075^{+1045}_{-1033}\text{ MeV}$, and $K_{\text{sym},0}=-156^{+97}_{-81}\text{ MeV}$ ($68\%$ standard deviations). This results in $90\%$ confidence intervals of 69 MeV $\leq K_{0}\leq$ 352 MeV, 1371 MeV $\leq M_{0}\leq$ 4808 MeV, and -285 MeV $\leq K_{\mathrm{sym,0}}\leq$ 7 MeV. We tabulate these values for both priors on $L_{0}$ in Table 1 for comparison to the simple method described in Sec. V.1. These constraints on the nuclear parameters are comparable to, yet smaller than that found in Sec. V.1, although are much more accurate because the covariances between $\tilde{\Lambda}$ and such nuclear parameters were properly taken into account, as well as considering the true probability distribution on $\tilde{\Lambda}$ from GW170817 as derived by the LIGO Collaboration.

How much does the addition of PEs affect the bounds on $K_{0}$, $M_{0}$, and $K_{\text{sym},0}$? To address this, we repeat our analysis without including these additional EoSs (see Appendix B for more details). We find that the removal of such EoSs gives strong improvement in both correlations and nuclear constraints for low-order nuclear parameters $K_{0}$ and $M_{0}$, and the results are consistent with those in Ref. Malik et al. (2018). This further illuminates the need to study a wider variety of EoSs for use in universal relations to properly account for systematic errors.

VI Conclusion and Discussion

The recent GW observation GW170817 coupled with the IR/UV/optical counterpart placed upper and lower bounds on the mass-weighted average tidal deformability $\tilde{\Lambda}$. We take advantage of this by selecting a diverse set of NS EoSs encompassing non-relativistic Skyrme-type interactions, RMF interactions and phenomenological variation of nuclear parameter models in order to constrain the nuclear matter parameters which are vital to limiting physically valid EoSs. We first found that approximate universal relations exist between linear combinations of nuclear parameters and $\tilde{\Lambda}$ for all values of mass ratio $q$ allowed from GW170817. We next constructed 2-dimensional probability distributions between $\tilde{\Lambda}$ and such nuclear parameters, converted them into one-dimensional conditional probability distributions on $K_{\text{sym,0}}$ given observations of $\tilde{\Lambda}$, and finally combined them with a posterior probability distribution on $\tilde{\Lambda}$ from LIGO and integrated them over $\tilde{\Lambda}$ in order to obtain posterior distributions on the nuclear parameters. From these posterior distributions, we derived 90% confidence intervals on the incompressibility $K_{0}$, its slope $M_{0}$, and the curvature of symmetry energy $K_{\text{sym},0}$ at saturation density as 69 MeV $\leq K_{0}\leq$ 352 MeV, 1371 MeV $\leq M_{0}\leq$ 4808 MeV, and -285 MeV $\leq K_{\mathrm{sym,0}}\leq$ 7 MeV. The bounds on $M_{0}$ and $K_{\text{sym},0}$ are more conservative and safer to quote than those found in Malik et al. (2018). In addition, the constraints derived on $K_{\text{sym,0}}$ shows agreement with those in Refs. Margueron and Gulminelli (2018); Mondal et al. (2018). We also note that bounds on $K_{0}$ and $M_{0}$ are less reliable than those on $K_{\mathrm{sym,0}}$ due to smaller correlations in the universal relations.

The bounds derived in this paper are only valid for NSs and may not be valid for hybrid stars (HSs) with quark core and nuclear matter envelope. We discuss this point in more detail in Appendix C.

Future work on this subject includes investigation into combinations of nuclear parameters other than the linear ones studied here, to see if the correlations among such new combinations against $\tilde{\Lambda}$ improves. For example, one can consider “multiplicative” combinations of the form $K_{0}L_{0}^{\eta}$ with constant $\eta$, in a similar spirit to Sotani et al. (2014); Silva et al. (2016). Furthermore, one can study how the universal relations considered in this paper change as a function of the chirp mass, which may be useful for future binary NS merger events. We also plan to study how the bounds derived here on nuclear parameters will improve in the future by considering upgraded ground-based GW detectors, such as aLIGO with its design sensitivity aLI , A+ Ap_ , Voyager Ap_ , Einstein Telescope ET and Cosmic Explorer Ap_ , in particular by combining multiple events, and at what point systematic errors due to the EoS variation in the universal relations dominate statistical errors on $\tilde{\Lambda}$. Work along these directions is currently in progress Carson et al. (2019).

Acknowledgments

We thank David Nichols for his illuminating advice on conditional probability distributions. K.Y. acknowledges support from NSF Award PHY-1806776. K.Y. would like to also acknowledge networking support by the COST Action GWverse CA16104. A.W.S. was supported by NSF grant PHY 1554876 and by the U.S. DOE Office of Nuclear Physics.

Appendix A $\tilde{\Lambda}$ versus $\Lambda_{1.4}$

Malik et al. Malik et al. (2018) first studied correlations between nuclear and tidal parameters for individual NSs. Given that the tidal parameter measured from GW observations is $\tilde{\Lambda}$, corresponding to the mass-weighted average of two tidal parameters in a binary, the authors of Ref. Malik et al. (2018) assumed the masses of the two NSs in GW170817 to be $m_{1}=1.4M_{\odot}$ and $m_{2}=1.33M_{\odot}$. Next they studied correlations between $\tilde{\Lambda}$ in such a binary and $\Lambda_{1.4}$, representing the tidal deformability for an individual NS with a mass of $1.4M_{\odot}$.

The above assumption can be dangerous because the individual mass measurements of GW170817 are not very accurate. Although the chirp mass has been measured with great accuracy as $\mathcal{M}=(m_{1}m_{2})^{3/5}(m_{1}+m_{2})^{-1/5}=1.188^{+0.004}_{-0.002}\text{ M}_{\odot}$, the mass ratio varies as $q=m_{2}/m_{1}\in[0.73,1.00]$ Abbott et al. (2017b).

The top panel of Fig. 11 presents the $\tilde{\Lambda}$–$\Lambda_{1.4}$ correlation for various $q$ within the above range with the chirp mass fixed to $\mathcal{M}=1.188M_{\odot}$, while the bottom panel shows the absolute fractional difference from the linear fit. Observe that a strong correlation exists between $\tilde{\Lambda}$ and $\Lambda_{1.4}$ for any $q$. The maximum fractional error for this case is $\sim 5$%, with a correlation coefficient of $C=0.998$. On the other hand, once we include the hybrid EoSs discussed in more detail in Appendix C, one clearly sees a large deviation from the correlation with other EoSs, with the fractional difference reaching up to 60%.

The behavior in Fig. 11 can be understood from Fig. 12, where we show $\tilde{\Lambda}$ against $q$ with $\mathcal{M}$ fixed to the measured value for GW170817. If we do not consider hybrid EoSs, $\tilde{\Lambda}$ is insensitive to $q$ Radice et al. (2018); Burgio et al. (2018), which is the origin of the strong correlation in the $\tilde{\Lambda}$–$\Lambda_{1.4}$ relation. On the other hand, for hybrid EoSs considered here, GW170817 can be either HS/HS or HS/NS when the mass ratio is close to unity555We note that hybrid EoSs considered in Paschalidis et al. (2018) admit either NS/NS or HS/NS for GW170817.. Thus, one finds a significant drop in $\tilde{\Lambda}$ as one increases $q$ Burgio et al. (2018), which changes the $\tilde{\Lambda}$–$\Lambda_{1.4}$ relation drastically.

Appendix B Repeated Analysis without PEs

In this appendix, we study the effect of PEs on nuclear parameter bounds by re-analyzing them without including such EoSs. This way, we can directly compare our results with those in Ref. Malik et al. (2018) which did not include these additional EoSs. Figure 13 once again presents correlations between $\tilde{\Lambda}$ and linear combinations of nuclear parameters as a function of mass ratio. Here, for comparison purposes we choose $\alpha=1.10$, $\beta=15.62$, and $\gamma=2.81$ such that correlations become maximum, as was done in Ref. Malik et al. (2018). Observe that correlations with $\tilde{\Lambda}$ remain almost constant throughout the entire region of allowable mass ratios. In addition, note how correlations for linear combination involving $K_{0}$ and $M_{0}$ are increased by up to $55\%$ from Fig. 5 which includes PEs, while linear combinations with higher order nuclear parameter $K_{\text{sym},0}$ interestingly shows a small decrease in correlation, yet remains comparable. This is revealing of the flexible nature of the $K_{\text{sym},0}$ nuclear parameter.

We now derive constraints on nuclear parameters without PEs. Following the procedure outlined in Sec. V, new bounds on $K_{0}$, $M_{0}$, and $K_{\text{sym},0}$ are calculated for a central mass ratio of $q=0.87$, and summarized in Table 2. Comparing this with Table 1, one sees that the additional PEs significantly weaken estimated constraints for low order nuclear parameters $K_{0}$ and $M_{0}$, and interestingly, improve them for high-order nuclear parameter $K_{\mathrm{sym,0}}$. Here we find results somewhat agreeable to what was found in Malik et al Malik et al. (2018), however enlarged due to the addition of EoS variation uncertainties.

Appendix C Hybrid Quark-hadron Stars

In this appendix, we investigate the use of an additional valid class of EoS: hybrid quark-hadron stars based on Ref. Paschalidis et al. (2018). Here, the low-density nucleonic matter region of PEs transition into a high-density quark matter phase in a given transitional energy density region $\epsilon_{1}\leq\epsilon\leq\epsilon_{2}$. For our purposes, we consider Set I quark matter EoSs, where the pressure following transition is given by Alford and Sedrakian (2017) (see also Montana et al. (2019); Seidov (1971); Zdunik and Haensel (2013); Alford et al. (2013)):

[TABLE]

with $c_{s}$ being the constant speed of sound in the quark matter, $\epsilon_{1}$ and $\epsilon_{2}$ characterizing the energy density “jump” $\epsilon_{2}-\epsilon_{1}\equiv\epsilon_{1}j$, and $P_{\text{tr}}$ representing the transition pressure, such that the low density hadronic matter’s energy density equals $\epsilon_{1}$. In this paper, we adopt the ACS-II parameterization in Paschalidis et al. (2018) as $P_{\text{tr}}=1.7\times 10^{35}$dyn/cm2, $\epsilon_{2}=8.34\times 10^{14}$g/cm3 and $c_{s}^{2}=0.8$ with $j=0.8$ or 1.

As we show in Fig. 12, strong phase transitions in the star admit a secondary stable HS configuration (denoted HS/HS). HSs evaluate to a reduction in tidal deformability $\tilde{\Lambda}$ from their NS-branch counterparts, thus altering universal relations accordingly. Here, we examine how this additional possibility of binary HSs and the choice of fiducial nuclear matter EoS impacts correlations between $\tilde{\Lambda}$ and nuclear parameters.

Figure 14 investigates this phenomena by choosing 3 different fiducial nuclear matter EoSs with soft ($\tilde{\Lambda}\approx 465$), intermediate ($\tilde{\Lambda}\approx 800$), and stiff ($\tilde{\Lambda}\approx 1045$) representative values of tidal deformability for $\mathcal{M}=1.188M_{\odot}$ and $q=0.995$. Next, HS EoSs are formulated, and new universal relations are derived - including both stellar configurations at high values of $q$, as can clearly be seen by the dashed vertical line in Fig. 12. Observe how the choice of fiducial nuclear matter EoS impacts the universal relations differently depending on which combination of nuclear parameters is used. For example, use of the stiff fiducial EoS compared to the intermediate one results in a small decrease in correlation for $K_{0}+\alpha L_{0}$, a negligible decrease for $M_{0}+\beta L_{0}$, and a large decrease for $K_{\text{sym},0}+\alpha L_{0}$. Alternatively, choice of the soft fiducial EoS results in medium decreases in correlation for $K_{0}+\alpha L_{0}$ and $M_{0}+\alpha L_{0}$, and an increase in correlation for $K_{\text{sym},0}+\alpha L_{0}$.

In conclusion, we find that the use of valid hybrid quark-hadron star EoSs in universal relations can influence universality in unexpected ways. Thus, the bounds derived in Table 1 are strictly valid only for NSs, and they are subject to change once one includes the possibility for HSs. Refer also to Ref. Montana et al. (2019) for a more detailed analysis of hybrid star EoSs in conjunction with GW170817.