Future Prospects for Constraining Nuclear Matter Parameters with Gravitational Waves

TL;DR

This paper explores how future gravitational wave observations from neutron star mergers can better constrain nuclear matter parameters by analyzing correlations between tidal deformabilities and nuclear physics, considering multiple events and systematic errors.

Contribution

It extends previous correlation analyses to varying chirp masses and proposes methods to improve nuclear parameter constraints using multiple detections and multidimensional correlations.

Findings

Future GW observations can improve nuclear parameter bounds by ~30%.

Systematic errors mainly stem from uncertainties in the equation of state.

Multidimensional correlations can reduce systematic errors.

Abstract

The gravitational wave emission from the merging binary neutron star system GW170817 arrived full of tidal information which can be used to probe the fundamental ultra-dense nuclear physics residing in these stars. In previous work, we used two-dimensional correlations between nuclear matter parameters and tidal deformabilities of neutron stars applying specifically to GW170817 to derive constraints on the former. Here, we extend this analysis by finding similar correlations for varying chirp masses, the dominant determining factor in the frequency evolution of the inspiral, such that one can apply the same method to future detections. We estimate how accurately one can measure nuclear parameters with future gravitational wave interferometers and show how such measurements can be improved by combining multiple events. We find that bounds on the nuclear parameters with future observations can improve from the current one with GW170817 only by $\sim 30\%$ due to the existence of systematic errors caused mainly by the remaining uncertainty in the equation of state near and just above the nuclear saturation density. We show that such systematic errors can be reduced by considering multidimensional correlations among nuclear parameters and tidal deformabilities with various neutron star masses.

Figures

Tables

Table 1. Table 1: Respective fitting functions for the covariance matrix 𝚺 𝚺 \bm{\Sigma} and the mean vector 𝝁 𝝁 \bm{\mu} in Eq. ( 16 ) necessary for the full reconstruction of the two-dimensional probability distributions between Λ ~ ~ Λ \tilde{\Lambda} and K sym , 0 subscript 𝐾 sym 0 K_{\text{sym},0} . Here, the values of Σ KK subscript Σ KK \Sigma_{\text{KK}} and μ K subscript 𝜇 K \mu_{\text{K}} correspond to the variance and mean of K sym , 0 subscript 𝐾 sym 0 K_{\text{sym},0} , which are independent of chirp mass, thus require no fitting function.

Parameter	Fitting Function
$\sqrt{{\Sigma }_{\Lambda \Lambda }}\left(ℳ\right)$ [MeV]	$\text{Exp}\left[6.287-11.86\log ℳ-0.9803{\log }^2ℳ\right]$
$\sqrt{{\Sigma }_{\Lambda \text{K}}}\left(ℳ\right)$ [MeV]	$\text{Exp}\left[5.073-5.477\log ℳ-4.103{\log }^2ℳ\right]$
$\sqrt{{\Sigma }_{\text{KK}}}\left(ℳ\right)$ [MeV]	$80.11$
${\mu }_{\Lambda }\left(ℳ\right)$ [MeV]	$\text{Exp}\left[7.394-12.27\log ℳ-6.399{\log }^2ℳ\right]$
${\mu }_{\text{K}}\left(ℳ\right)$ [MeV]	$-71.6$

Equations

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

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

$$\Phi (x^i)=$$

$$=$$

$$\displaystyle+\frac{1}{2}\mathcal{E}_{ij}x^{i}x^{j}+\mathcal{O}\Bigg{(}\frac{r^{3}}{R^{3}}\Bigg{)},$$

$$g_{\alpha \beta }=g_{\alpha \beta }^{(0)}+h_{\alpha \beta },$$

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

$$\mathcal{M}\equiv\Bigg{(}\frac{q^{3}}{1+q}\Bigg{)}^{1/5}m_{1}.$$

$$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}$$

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

$$\Sigma_{ab}=\frac{1}{N}\sum^{N}_{i=0}a_{i}b_{i}-\frac{1}{N^{2}}\Big{(}\sum^{N}_{i=0}a_{i}\Big{)}\Big{(}\sum^{N}_{i=0}b_{i}\Big{)}.$$

$$P_A(\tilde{\Lambda })=\frac{1}{\sqrt{2\pi {\sigma }_A^2}}{e}^{-(\tilde{\Lambda }-{\mu }_{\tilde{\Lambda }}{)}^2/2{\sigma }_A^2}$$

$${\theta }^a=(\ln \mathcal{A},{\varphi }_c,t_c,\ln {\mathcal{M}}_z,\ln \eta ,{\chi }_s,{\chi }_a,\tilde{\Lambda },\delta \tilde{\Lambda }),$$

$$\Delta\theta^{a}=\sqrt{\Big{(}\tilde{\Gamma}^{-1}\Big{)}^{aa}}.$$

$$\tilde{\Gamma}_{ab}\equiv\Big{(}\frac{\partial h}{\partial\theta^{a}}\Big{|}\frac{\partial h}{\partial\theta^{b}}\Big{)}+\frac{1}{\sigma^{2}_{\theta^{a}}}\delta_{ab},$$

$$(a|b)\equiv 2\underset{0}{\overset{\infty }{\int }}\frac{{\tilde{a}}^{\ast }\tilde{b}+{\tilde{b}}^{\ast }\tilde{a}}{S_n(f)}df.$$

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

$$\log y_i=a_i+b_i\log \mathcal{M}+c_i(\log \mathcal{M}{)}^2,$$

$$P(K_{\text{sym,0}}|\tilde{\Lambda}_{\text{obs}})\sim\mathcal{N}\Bigg{(}\mu_{K_{\text{sym,0}}}+\frac{\sigma_{K_{\text{sym,0}}}}{\sigma_{\tilde{\Lambda}_{\text{obs}}}}C(\tilde{\Lambda}_{\text{obs}}-\mu_{{\tilde{\Lambda}_{\text{obs}}}}),(1-C^{2})\sigma_{K_{\text{sym,0}}}^{2}\Bigg{)}.$$

$$P_A(K_{\text{sym},0})=\underset{-\infty }{\overset{\infty }{\int }}P(K_{\text{sym},0}|\tilde{\Lambda })P_A(\tilde{\Lambda })d\tilde{\Lambda },$$

$$\Lambda \approx {\Lambda }_{1.4}+{\Lambda }_{1.4}^{\prime }\left(1-\frac{m}{m_0}\right).$$

$$\sqrt{\mathbf{\Sigma }}=\left(\begin{array}{cc}193.6& 97.10\\ 97.10& 80.11\end{array}\right)[\mathrm{MeV}],\mathbf{\mu }=\left(\begin{array}{c}543.2\\ -71.164\end{array}\right)[\mathrm{MeV}],$$

$$P^{\prime \prime \prime }(K_{\text{sym},0})=P(K_{\text{sym},0},{\bar{\Lambda }}_x,{\bar{\Lambda }}_y,{\bar{\Lambda }}_{1.5}).$$

$$P_{\text{O2}}(\tilde{\Lambda })=\frac{1}{\sqrt{2\pi (172.5{)}^2}}{e}^{-(\tilde{\Lambda }-430.8{)}^2/2(172.5{)}^2}.$$

$$\mathbf{\Sigma }=\left(\begin{array}{cc}45610& 10410\\ 10410& 6418\end{array}\right),\mathbf{\mu }=\left(\begin{array}{c}606.7\\ -71.16\end{array}\right),$$

$$P(K_{\text{sym,0}}|{\tilde{\Lambda }}_{\text{obs}})=\frac{\text{Exp}\left[-\frac{{\left(K_{\text{sym,0}}-(0.228({\tilde{\Lambda }}_{\text{obs}}-606.7)-71.16)\right)}^2}{2(63.59{)}^2}\right]}{\sqrt{2\pi (63.59{)}^2}}.$$

$$P(K_{\text{sym},0})=\underset{-\infty }{\overset{\infty }{\int }}P(K_{\text{sym},0}|\tilde{\Lambda })P_{\text{O2}}(\tilde{\Lambda })d\tilde{\Lambda }.$$

$$N_A=\Delta {\tau }_0\underset{0}{\overset{z_h}{\int }}4\pi [a_0{r}_1(z){]}^2\mathcal{R}r(z)\frac{d\tau }{dz}dz.$$

$$a_0{r}_1(z)=\frac{1}{H_0}\underset{0}{\overset{z}{\int }}\frac{d{z}^{\prime }}{\sqrt{(1-{\Omega }_{\Lambda })(1+z^{\prime }{)}^3+{\Omega }_{\Lambda }}},$$

$$\frac{d\tau }{dz}=\frac{1}{H_0}\frac{1}{1+z}\frac{1}{\sqrt{(1-{\Omega }_{\Lambda })(1+z^{\prime }{)}^3+{\Omega }_{\Lambda }}},$$

Full text

I Introduction

Neutron stars (NSs) exist in one of the most extreme states of matter found in the universe. However, the determination of the equation of state (EoS) of ultra-dense matter found exclusively in such compact objects remains to be one of the largest unsolved mysteries in both nuclear physics and astrophysics to date. The nuclear matter EoS determines many important stellar properties, such as the mass, radius, and tidal deformability, and is vital to the further study of supranuclear matter. Independent measurements of certain macroscopic NS observables determined by the EoS, such as the mass and radius, can be used to constrain the EoS, as was indeed done in Refs. Guver and Ozel (2013); Ozel et al. (2010); Steiner et al. (2010); Lattimer and Steiner (2014); Ozel and Freire (2016) via x-ray observations of the mass-radius relationship. However, such measurements potentially suffer from large systematic errors due to uncertainties in the astrophysical modeling of x-ray bursts.

Recent observations of gravitational waves (GWs) from a merging binary NS system (GW170817 Abbott et al. (2017a)) have been used to probe the interior nuclear structure via imprinted tidal effects Abbott et al. (2019, 2018a); Paschalidis et al. (2018); Burgio et al. (2018); Malik et al. (2018); Landry and Essick (2019), which offers us a cleaner method of determining the nuclear matter EoS than electromagnetic wave observations of neutron stars. As the NSs lose energy via GW emission, they inspiral towards each other and become increasingly tidally deformed in response to the companion stars’ tidal field. This deformation is characterized by the tidal deformability Flanagan and Hinderer (2008) of the NS, and is strongly dependent on the underlying EoS. Further, the mass-weighted combination of such tidal deformabilities from each star is the leading tidal parameter in the gravitational waveform, which has been constrained by the LIGO and VIRGO Collaboration to a 90% credible bound of $70\leq\tilde{\Lambda}\leq 720$ Abbott et al. (2019); De et al. (2018). Such observations have also been mapped to the NS radius in Refs. Annala et al. (2018); Abbott et al. (2018a); Lim and Holt (2018); Bauswein et al. (2017); De et al. (2018); Most et al. (2018).

While all currently-proposed EoSs to date utilize various different approximations, one way to effectively study them is by measuring the nuclear parameters which parameterize the EoSs using a model-independent formalism. One such method for doing this111Piecewise polytropic constructions 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. (2018a) similarly parameterize nuclear matter EoSs in a model-independent way. See also Landry and Essick (2019); Kumar and Landry (2019) for related works on piecewise unified EoSs. is to Taylor expand the energy per nucleon of asymmetric nuclear matter about the saturation density of symmetric matter Myers and Swiatecki (1969); Vidaña et al. (2009). The resulting coefficients are known as the “nuclear matter parameters” and consist of the following: the slope of the symmetry energy $L_{0}$; the nuclear incompressibility $K_{0}$; the slope of the incompressibility $M_{0}$; the curvature of the symmetry energy $K_{\text{sym},0}$; and higher orders, each evaluated at the nuclear saturation density.

Previous important analyses by Alam et al. Alam et al. (2016) found approximately universal relations between the NS radius at a given mass and the nuclear parameters mentioned above (similar work can be found in Refs. Sotani et al. (2014); Silva et al. (2016)). Further, Malik et al. Malik et al. (2018) found that certain linear combinations of nuclear parameters (such as $K_{0}+\alpha L_{0}$ with $\alpha$ chosen to give maximal correlation) gave way to heightened correlations with the individual tidal deformabilities evaluated at a given mass. By assuming individual masses for GW170817 and taking the approximate universal relations to be exact, the authors utilized prior constraints on the tidal deformability from GW170817 Abbott et al. (2017b); Radice et al. (2018) and $L_{0}$ Abbott et al. (2019); Oertel et al. (2017); Lattimer and Steiner (2014) to derive new constraints on the nuclear parameters.

In Ref. Carson et al. (2019a), we expanded upon the previous work of Malik et al. Malik et al. (2018) and found improved constraints on $K_{0}$, $M_{0}$, and $K_{\text{sym},0}$. We made various improvements upon the original analysis, such as (i) using an expanded set of EoSs in order to fully take into account the EoS-variation systematic errors; (ii) considering the correlations with the mass-weighted tidal deformability $\tilde{\Lambda}$ (rather than $\Lambda_{1.4}$, the tidal deformability at $1.4\text{ M}_{\odot}$) which was directly measured by GW observations; (iii) fully considering the EoS variation scatter uncertainty in the estimation of constraints; and (iv) adopting the full posterior distribution on $\tilde{\Lambda}$ as measured from GW170817. Further, we found that by using our method of computing the posterior distribution on nuclear parameters, high degrees of correlation between $\tilde{\Lambda}$ and such nuclear parameters was not necessary, as any covariances between the two was taken into account by the multi-variate Gaussian probability distribution between them. On the other hand, the authors of Ref. Malik et al. (2018) absolutely required high correlations with nuclear parameters in order to assume that the relationship between the tidal deformability and the nuclear parameters lie exactly on the best-fit line between the two. The resulting 90% confidence interval on the curvature of symmetry energy was found to be $-259\text{ MeV}\leq K_{\text{sym},0}\leq 32\text{ MeV}$, which was more conservative than that derived in Ref. Malik et al. (2018).

Here, we extend the work of Ref. Carson et al. (2019a) into the future of GW astronomy. While every future merger event will be composed of NSs with varying individual masses which are difficult to measure, we can categorize them by the chirp mass $\mathcal{M}$, which is the dominant driving factor in the frequency evolution of the inspiral event given by a certain combination of individual NS masses. In this investigation, we repeat the analysis of Ref. Carson et al. (2019a) as a function of chirp mass, applicable to any future event. Further, we restrict the set of EoSs to those that obey the nuclear parameter correlations of Ref. Tews et al. (2017), and consider the implications of observations using future GW interferometers: Advanced LIGO (aLIGO) aLI , LIGO A+ (A+) Ap_ , Voyager Ap_ , Cosmic Explorer (CE) Ap_ , and Einstein Telescope (ET) ET . We will consider not only the increased sensitivities from current detectors but also the combined uncertainties from multiple-event detections (relevant for future detectors with expanded horizon volumes).

I.1 Executive summary

Here we summarize our results for busy readers. In the current analysis, we attempt to find constraints on the nuclear matter parameters as a function of the binary systems’ chirp mass. We begin by finding the correlations between the mass-weighted tidal deformability $\tilde{\Lambda}$ and various nuclear parameters, and combinations thereof. Similar to our results in Ref. Carson et al. (2019a), we find that the low-order nuclear parameters $K_{0}$ and $M_{0}$ observe small correlations with $\tilde{\Lambda}$, while the correlations for higher-order parameter $K_{\text{sym},0}$ remains high at $\sim 80\%$. For this reason, we consider constraints on the curvature $K_{\text{sym},0}$ of the symmetry energy which is one of the most uncertain parts of the EoS of dense nucleonic matter Li et al. (2019).

In this analysis, we compute the posterior probability distribution for the curvature of symmetry energy $K_{\text{sym},0}$, for 22 different values of chirp mass between $0.94\text{ M}_{\odot}$ and $1.6\text{ M}_{\odot}$. For each value of chirp mass $\mathcal{M}_{i}$ considered, we compute the single-event $\tilde{\Lambda}$ uncertainties using Fisher analysis techniques. Approximated as a Gaussian prior, the uncertainty in $\tilde{\Lambda}$ may be used to estimate the posterior probability distribution on $K_{\text{sym},0}$, by multiplication with the one-dimensional conditional probability distribution on $K_{\text{sym},0}$ given an observation of $\tilde{\Lambda}$, and then integrating over all $\tilde{\Lambda}$ values. The process is then repeated for each value of chirp mass $\mathcal{M}_{i}$, resulting in a relationship between the uncertainties in $K_{\text{sym},0}$ and the chirp mass $\mathcal{M}$.

The corresponding one-sided 90% confidence intervals on $K_{\text{sym},0}$ for single-event detections, along with the calculated systematic errors, are plotted in Fig. 1 as a function of chirp mass for 6 different GW interferometers. In this figure, we observe that as the detector sensitivity is increased, the statistical errors become subdominant rather quickly. This in turn forces the overall errors to approach the systematic error “wall” at $\sim 104$ MeV, caused partially by uncertainties in the EoS at low-density which are less sensitive to neutron star tidal deformabilities. For this reason, the curves corresponding to the 3rd-generation detectors CE and ET become indistinguishable from the systematic errors – indicating the necessity to reduce such errors in order for the further constraint of $K_{\text{sym},0}$ to become possible.

Following this, we offer a method to further decrease the statistical errors in the measurement of $K_{\text{sym},0}$. This is accomplished by repeating the same analysis for the fixed chirp mass of $1.188\text{ M}_{\odot}$ with the first coefficient, $\Lambda_{1.4}$, of the Taylor expanded tidal deformability $\Lambda\approx\Lambda_{1.4}+\Lambda_{1.4}^{\prime}(1-\frac{m}{m_{0}})$ about $m_{0}=1.4\text{ M}_{\odot}$222$m_{0}$ will remain fixed for the remainder of the analysis, with the exception of Sec. V, where we consider the effect of variations in $m_{0}$., rather than the mass-weighted tidal deformability $\tilde{\Lambda}$. $\Lambda_{1.4}$ (or the tidal deformability at $1.4\text{ M}_{\odot}$) is mass-independent, and thus, it is identical for all future GW events. This way, we can combine the uncertainties for multiple detected events when it becomes applicable for future detectors. As was observed in Fig. 1, we found that the uncertainties in $K_{\text{sym},0}$ became dominated by systematics for the single-event analyses on Voyager-era detectors and beyond. By combining GW170817-like events detected on aLIGO and A+, we find that one can further reduce the statistical errors in $K_{\text{sym},0}$ such that the errors similarly become dominated by systematics.

Finally, we investigate the reduction of systematic errors by adding information about the tidal deformability at various different masses. We begin by generating a four-dimensional Gaussian probability distribution $P(K_{\text{sym},0},\Lambda_{m_{x}},\Lambda_{m_{y}},\Lambda_{m_{z}})$. The systematic error on $K_{\text{sym}}$ is obtained by first evaluating $P(K_{\text{sym},0},\Lambda_{m_{x}},\Lambda_{m_{y}},\Lambda_{m_{z}})$ at the fiducial values of $\Lambda_{m_{x}}$, $\Lambda_{m_{y}}$, and $\Lambda_{m_{z}}$ and estimating the 90% confidence interval. Figure 2 presents the systematic errors when $m_{z}$ is fixed to be $1.5\text{ M}_{\odot}$, and ($m_{x},m_{y}$) are varied between $1.0\text{ M}_{\odot}$ and $2.0\text{ M}_{\odot}$. Observe how the resulting systematic uncertainties are reduced to $\sim 74$ MeV for certain combinations of $m_{x}$ and $m_{y}$.

The remainder of this paper is organized as follows. We begin in Sec. II with a review on the tidal deformability of neutron stars, as well as the nuclear matter EoS and its constituent nuclear parameters. We follow this up in Sec. III with a study on the correlations between various nuclear parameters and the mass-weighted tidal deformability. In Sec. IV, we analyze the measurement accuracy on such nuclear parameters as a function of the chirp mass, applicable to future binary NS mergers. In Sec. V, we discuss how one can further reduce systematic errors by considering correlations among nuclear parameters and multiple tidal deformabilities at different NS masses. We conclude in Sec. VI with a discussion on the implications of our results, as well as an outline into possible avenues of future work. We have adopted geometric units such that $G=c=1$ throughout.

II Background and theory

In this section we begin with a review on the NS tidal deformability in Sec. II.1, followed up by a review on the NS equation of state and its constituent nuclear parameters in Secs. II.2 and II.3.

II.1 Neutron star tidal deformability

Here we offer a brief overview on how one can extract information on the internal structure of a NS by way of the GW observations of binary NS merger events, such as GW170817. In the presence of a neighboring tidal field $\mathcal{E}_{ij}$ (such as a NS in a binary system with a compact companion), NSs will acquire a quadrupole moment $Q_{ij}$ characterized as the linear response to $\mathcal{E}_{ij}$:

[TABLE]

with tidal deformability $\lambda$ Flanagan and Hinderer (2008); Hinderer (2008); Damour and Nagar (2009); Binnington and Poisson (2009); Yagi and Yunes (2013). The tidal deformability characterizes the NSs corresponding deformation from sphericity, and can be made unitless by the following normalization:

[TABLE]

with stellar mass $M$.

Following Refs. Hinderer (2008); Damour and Nagar (2009); Yagi and Yunes (2013); Carson et al. (2019a), the dimensionless tidal deformability $\Lambda$ can be computed by isolating different asymptotic limits of the gravitational potential in the buffer zone $R\ll r\ll\mathcal{L}$ given by

[TABLE]

where $r=|x^{i}|$, $\mathcal{L}$ is the length scale of the companion-induced curvature, and $R$ is the stellar radius. Here, $g_{tt}$ corresponds to the $tt$-component of the full spacetime metric:

[TABLE]

constructed via a non-spinning, spherically-symmetric background solution $g_{\alpha\beta}^{(0)}$ perturbed by the tidal deformation with metric components $h_{\alpha\beta}$. The perturbed Einstein equations may then be solved in the interior of the NS and matched to the exterior solution at the surface of the star. Further, the radius $R$ may then be determined from the condition $p(R)=0$.

In this investigation, we consider the scenario of two NSs orbiting each other in a binary system, like GW170817. In this case, each NS individually obtains quadrupole moments from the neighboring tidal field, resulting in two highly-correlated tidal deformabilities $\Lambda_{1}$ and $\Lambda_{2}$. Due to these correlations, individual tidal deformabilities are very difficult to extract from GW observations. Typically, it is useful to reparameterize the waveform via independent linear combinations of $\Lambda_{1}$ and $\Lambda_{2}$ which enter the gravitational waveform at 5th post-Newtonian (PN) and 6PN orders333$n$PN order corrections enter the gravitational waveform at relative powers of $(v/c)^{2n}$. respectively. The dominant tidal effect in the resulting waveform is known as the mass-weighted tidal deformability, and is given by Flanagan and Hinderer (2008)

[TABLE]

with mass ratio $q\equiv m_{2}/m_{1}$ ($m_{1}\geq m_{2}$). Here we also define the chirp mass of the binary system, which is the primary controlling factor of the merger inspiral defined by

[TABLE]

Similarly to $\tilde{\Lambda}$, this quantity can be measured with much higher accuracy than either of the individual masses $m_{1}$, $m_{2}$, or the mass ratio $q$. For this reason, we consider the binary chirp mass $\mathcal{M}$ to be the dominant dependent variable in this analysis, cataloging our various results as a function of $\mathcal{M}$ for any future GW event.

II.2 The nuclear matter parameters

While the NS EoS is not currently known, there are many methods one can use to restrict it using various observations. This is because the structure of a NS and many of its observables such as mass, radius, tidal deformability, etc. rely strongly on the underlying EoS of nuclear matter. For example, GW observations may help constrain the EoS in the pressure-density plane Abbott et al. (2018b). In this paper, we show how GW detections can aid in the constraint of various characteristics of the EoS, known as the nuclear matter parameters Vidaña et al. (2009).

As originally considered in Ref. Myers and Swiatecki (1969) and followed up in Refs. Vidaña et al. (2009); Alam et al. (2016); Malik et al. (2018); Carson et al. (2019a), we offer a generic method to parameterize NS EoSs. This is done by first Taylor expanding the energy per nucleon $e$ of asymmetric nuclear matter about $\delta=0$ (symmetric nuclear matter), where $\delta\equiv(n_{n}-n_{p})/n$ is the isospin symmetry parameter for nuclear matter with $n_{n}$ neutron density, $n_{p}$ proton density, and total density $n=n_{n}+n_{p}$:

[TABLE]

Here $e(n,0)$ and $S_{2}(n)$ are the symmetric and second-order asymmetric nuclear matter energy per nucleon, respectively. Such energies are then further Taylor expanded about the nuclear saturation density $n_{0}\approx 2.3\times 10^{14}\text{ g/cm}^{3}$ as

[TABLE]

with $y\equiv(n-n_{0})/3n_{0}$. The above coefficients, all evaluated at the nuclear saturation density, determine the NS EoS and are referred to in the literature Myers and Swiatecki (1969); Vidaña et al. (2009) as the energy per particle $e_{0}$; incompressibility coefficient $K_{0}$; third derivative of symmetric matter $Q_{0}$; the slope of the incompressibility $M_{0}\equiv Q_{0}+12K_{0}$ Alam et al. (2016); Malik et al. (2018); symmetry energy $J_{0}$; its slope $L_{0}$; and its curvature $K_{\mathrm{sym,0}}$.

In this paper, we expand upon previous works Alam et al. (2016); Malik et al. (2018); Carson et al. (2019a) and investigate the correlations between the mass-weighted tidal deformability $\tilde{\Lambda}$ and various nuclear parameters $L_{0}$, $K_{0}$, $M_{0}$, and $K_{\text{sym},0}$ in order to derive constraints on such parameters. Specifically, we focus on the curvature of the symmetry energy $K_{\text{sym},0}$, shown to be one of the most uncertain features of the nuclear matter EoS, especially at supranuclear densities found primarily in NSs Li et al. (2019).

II.3 The supranuclear equation of state

We now explain which EoSs we use in this paper. In Ref. Carson et al. (2019a), we showed the importance of considering a wide range of physically valid EoSs when computing constraints on nuclear parameters, in order to more properly take into account the systematic errors. In the current analysis, we employ a restricted set of the same EoSs as was used previously, taking into account the observed correlations between nuclear parameters $J_{0}$ and $L_{0}$. Starting with the 121 nuclear EoS models found in Ref. Carson et al. (2019a), we further remove 63 EoS models which do not comply with the allowed regions shown in Fig. 8 of Tews et al. Tews et al. (2017). Here, they combined an exclusion region $J_{0}(L_{0})$ with the “accepted” 95.4% correlation confidence bands between $J_{0}$ and $L_{0}$. See App. A for the impact of restricted EoSs on correlations between $K_{\mathrm{sym,0}}$ and $\tilde{\Lambda}$.

Taking the shared region between the above two exclusions results in 58 different nuclear EoS models, which can be classified into 3 distinct classes: 13 non-relativistic “Skyrme-type” EoSs, 5 relativistic-mean-field (RMF) EoSs, and 40 EoS models developed with a phenomenological variation method (PEs). The Skyrme-type models used consist of: SKa, Sly230a Chabanat et al. (1997), Sly2, Sly9 Chabanat (1995), Sly4 Chabanat et al. (1998), SkOp Reinhard (1999), SK255, SK272 Agrawal et al. (2003), BSK20, BSK21 Goriely et al. (2010), BSK22, BSK24, BSK26 Goriely et al. (2013). Further, the RMF models used are: BSR2, BSR6 Dhiman et al. (2007); Agrawal (2010), NL3$\omega\rho$ Carriere et al. (2003), DD2 Typel et al. (2010), DDH$\delta$ Gaitanos et al. (2004). All 18 of the above EoSs originate from the minimal set of EoSs used in Refs. Alam et al. (2016); Malik et al. (2018), now restricted by nuclear matter correlations. Lastly, following Ref. Read et al. (2009), the high-density core-region of the EoSs used in our analysis are matched to the low-density crust-region of the SLy EoS model Douchin and Haensel (2001) at about half of the nuclear saturation density, $\rho_{\text{stitch}}\approx 1.3\times 10^{14}\text{ g/cm}^{3}$.

One last class of EoSs indirectly used in our analysis can be found in Ref. Carson et al. (2019b), which we call “LVC constrained” EoSs in this paper. By sampling the full physical EoS parameter space, the LIGO and Virgo Collaboration Abbott et al. (2018b); Carney et al. (2018) derived a marginalized 90% posterior region on the NS pressure as a function of mass density (EoS) from GW170817, as seen in Fig. 2 of Abbott et al. (2018b). By randomly sampling the EoS posteriors from this analysis, a set of 100 “constrained” EoSs were obtained, restricted by the GW observation of GW170817. While we do not directly utilize these 100 EoSs in the current analysis, we use them to estimate the mean value of the mass-weighted tidal deformability $\tilde{\Lambda}$ in Sec. IV.1, seen by Fig. 6.

III Correlations between $\tilde{\Lambda}$ and nuclear parameters

We begin the current analysis by studying the correlations between the mass-weighted tidal deformability $\tilde{\Lambda}$ and various nuclear parameters $K_{0}$, $M_{0}$, and $K_{\text{sym},0}$ as was done in Ref. Carson et al. (2019a). In the previous analysis, we studied the correlations and constraints as a function of the mass ratio $q\equiv m_{1}/m_{2}$ for the fixed chirp mass of $1.188\text{ M}_{\odot}$; corresponding to GW170817. We here supplement this investigation by considering the correlations and constraints as a function of varying chirp mass at a fixed mass ratio, applicable to any number of future GW observations. The left panel of Fig. 3 justifies the use of a fixed mass ratio by presenting $\tilde{\Lambda}$ as a function of $q$ for the various EoS models used in this analysis. Observe how $\tilde{\Lambda}$ is insensitive to the choice of $q$. Such feature is absent in the right panel of Fig. 3, where $\tilde{\Lambda}$ is plotted as a function of chirp mass for fixed $q$. Thus, for the remainder of this analysis we fix the mass ratio to be $q=0.90$, corresponding to the center of $0.80\leq q\leq 1.00$ derived in Ref. Coughlin et al. (2019) for GW170817.

We measure the amount of correlation between two observables $x$ and $y$ via

[TABLE]

with covariances $\Sigma_{ab}$ $(a=(x,y),b=(x,y))$ given by

[TABLE]

Here $N$ represents the number of data points. A correlation of $\mathcal{C}=1$ represents perfect correlation between observables $x$ and $y$, while $\mathcal{C}=0$ corresponds to no correlation.

In Ref. Carson et al. (2019a), we studied the universal relations between $\tilde{\Lambda}$ and various nuclear parameters $K_{0}$, $M_{0}$, $K_{\text{sym},0}$ in order to derive constraints on such parameters. Additionally, it was shown in Refs. Alam et al. (2016); Malik et al. (2018); Carson et al. (2019a) that certain linear combinations of nuclear parameters, specifically $K_{0}+\alpha L_{0}$, $M_{0}+\beta L_{0}$, and $K_{\text{sym},0}+\gamma L_{0}$, exhibit heightened correlations, allowing one to derive more accurate constraints on the individual nuclear parameters. However, it was found in Ref. Carson et al. (2019a) that this came at the expense of additional sources of uncertainty which, if properly accounted for, enlarges the resulting constraints on the nuclear parameters $K_{0}$, $M_{0}$, and $K_{\text{sym},0}$. It was also found that the single nuclear parameters as well as the linear combinations involving $K_{0}$ and $M_{0}$ observed poor correlations of $\mathcal{C}\lesssim 0.50$; indicating somewhat unreliable constraints on the nuclear parameters.

Figure 4 similarly shows the above correlations as a function of chirp mass for a fixed mass ratio of $q=0.90$. Observe how, similar to what was found in Ref. Carson et al. (2019a), the correlations for $K_{0}$, $M_{0}$, $K_{0}+\alpha L_{0}$, and $M_{0}+\beta L_{0}$ are exceedingly poor for all values of chirp mass. $K_{\text{sym},0}$ on the other hand, remains highly correlated with $\tilde{\Lambda}$ across the entire range of $\mathcal{M}$. We also observe how correlations are not improved by much when considering linear combinations between $K_{\text{sym},0}$ and $L_{0}$.

Could other combinations of nuclear parameters give stronger correlations? To address this question, we further explore new combinations of nuclear parameters in App. B. In particular, we consider the “multiplicative” combinations of $K_{0}L_{0}^{\eta}$, $M_{0}L_{0}^{\nu}$, and $K_{\text{sym},0}L_{0}^{\mu}$ that is motivated from Refs. Sotani et al. (2014); Silva et al. (2016). We found that such new combinations do not offer any advantages in terms of correlations and constraints.

For the above reasons, we consider only the curvature of the symmetry energy $K_{\text{sym},0}$ for the remainder of the investigation, without combinations with other parameters which would otherwise introduce additional uncertainties in the computation of constraints.

IV Nuclear parameter constraints with future GW observations

Now that we have identified the high-correlation behavior of $K_{\text{sym},0}$, we proceed to compute projected bounds on the curvature of the symmetry energy as a function of chirp mass that is applicable to any future event. Additionally we offer the same analysis repeated for 5 anticipated future detector sensitivities $S_{n}(f)$ for detectors O2 aLI , aLIGO aLI , A+ Ap_ , Voyager Ap_ , ET ET and CE Ap_ (see Fig. 5), which would allow one to compute the corresponding posterior distribution on $K_{\text{sym},0}$ given an events’ chirp mass $\mathcal{M}$.

IV.1 Single events

Previously in Ref. Carson et al. (2019a), a posterior distribution on $\tilde{\Lambda}$ as derived from GW170817, was utilized in order to compute posterior distributions on the nuclear parameters. In this analysis of future observations however, no such distribution is available. To remedy this, we approximate the effective “future” posterior distribution on $\tilde{\Lambda}$ as a Gaussian probability distribution given by

[TABLE]

for detector $A$. Here, $\mu_{\tilde{\Lambda}}=\mu_{\tilde{\Lambda}}(\mathcal{M})$ is computed from the mean value of the “LVC constrained” EoSs Carson et al. (2019b) described in Sec. II.3 for each value of chirp mass, as shown by Fig. 6. Further, $\sigma_{A}$ is approximated via simple Fisher analyses (described below), which estimates the measurement accuracy on $\tilde{\Lambda}$ under the assumption of detector sensitivity $A$. Figure 7 presents $\sigma_{A}$ for all 6 detectors.

The Fisher analysis method Cutler and Flanagan (1994); Berti et al. (2005); Poisson and Will (1995) estimates the accuracy with which one can extract best-fit parameters $\theta^{a}$ given prior probability distributions $\sigma_{\theta^{a}}$ and a template waveform $h$. For this analysis, we consider the sky-averaged “IMRPhenomD” gravitational waveform template $h$ Khan et al. (2016); Husa et al. (2016) for point particles, modified by the 5PN and 6PN tidal corrections given in Ref. Wade et al. (2014). We utilize a template parameter vector $\theta^{a}$ consisting of

[TABLE]

where $\mathcal{A}=\frac{\mathcal{M}_{z}^{5/6}}{\sqrt{30}\pi^{2/3}D_{L}}$ is a normalized amplitude factor, $D_{L}$ is the luminosity distance to the event, $\mathcal{M}_{z}=(1+z)\mathcal{M}$ is the redshifted chirp mass, $\eta=m_{1}m_{2}/(m_{1}+m_{2})^{2}$ is the symmetric mass ratio, $\chi_{s,a}=\frac{1}{2}(\chi_{1}\pm\chi_{2})$ are the symmetric and anti-symmetric combinations of individual spins $\chi_{1,2}$, and $\delta\tilde{\Lambda}$ is a tidal parameter entering first at 6PN order. In this investigation, we utilize fiducial parameter values of $\phi_{c}=0$, $t_{c}=0$, $\chi_{s}=\chi_{a}=0$, $\delta\tilde{\Lambda}=0$, $\mathcal{M}$ and $\eta$ are chosen from the current chirp mass iteration with a mass ratio of $q=0.9$, and $\tilde{\Lambda}$ is computed from the mean value of the “LVC constrained” EoSs. Additionally, we impose Gaussian spin priors of $|\chi_{s,a}|<1$, and tidal priors of $0<\tilde{\Lambda}<3000$, and $|\delta\tilde{\Lambda}|<500$ Wade et al. (2014).

To obtain a posterior distribution, one needs to know both the likelihood and prior distributions. Assuming the prior distributions on template parameters $\theta^{a}$ are Gaussian444A more comprehensive Bayesian analysis assumes the more valid choice of uniform prior distributions., the resulting posterior distributions are also Gaussian with root-mean-squares of

[TABLE]

The Fisher matrix $\tilde{\Gamma}_{ab}$ is defined as

[TABLE]

where $\sigma_{\theta^{a}}$ are the prior root-mean-square estimates of parameters $\theta^{a}$, and the inner product $(a|b)$ is given by

[TABLE]

Now that all of the tools are prepared, we next compute the posterior distributions on $K_{\text{sym},0}$ using the Gaussian prior distributions on $\tilde{\Lambda}$ computed above as a function of chirp mass, for future detectors. Following the process used in Ref. Carson et al. (2019a), this is accomplished by first generating a two-dimensional Gaussian probability distribution between $K_{\text{sym},0}$ and $\tilde{\Lambda}$, taking into account the covariances between the two as

[TABLE]

Here $\bm{x}$ and $\bm{\mu}$ are the 2D vectors containing $(\tilde{\Lambda},K_{\text{sym},0})$ and their means respectively, and $\bm{\Sigma}$ is the covariance matrix with elements given by Eq. (10).

Let us now offer readers the means to fully reproduce the results of the above analysis for any future event by constructing a fit for $\bm{\mu}$ and $\bm{\Sigma}$ in terms of chirp mass $\mathcal{M}$. Based on the relations between the former and the latter as shown in Fig. 8, we create a fit in a logarithmic power expansion as

[TABLE]

with $y_{i}$ being the various parameters $\sqrt{\Sigma_{ab}}$ and $\mu_{a}$, and fitting coefficients $a_{i}$, $b_{i}$, and $c_{i}$ which are summarized in Table 1. Observe how well the fit agrees with the numerical data in Fig. 8.

Constraints on $K_{\text{sym},0}$ are extracted by first computing the conditional probability distributions on $K_{\text{sym,0}}$ given a tidal deformability observation of $\tilde{\Lambda}_{\text{obs}}$. By following Ref. Jensen (2007), we can generate the one-dimensional conditional probability distribution on $K_{\text{sym,0}}$ by taking

[TABLE]

Above, $\mathcal{N}(\mu,\sigma^{2})$ is the normal distribution with mean and variance $\mu$ and $\sigma^{2}$, while $\mu_{\mbox{\tiny A}}$ and $\sigma_{\mbox{\tiny A}}^{2}$ are the mean and variances of $K_{\text{sym,0}}$ and $\tilde{\Lambda}_{\text{obs}}$. Finally, we can combine the one-dimensional conditional probability distribution function of Eq. (18) with the one-dimensional prior distribution on $\tilde{\Lambda}$ of Eq. (11). Marginalizing over $\tilde{\Lambda}$ results in a posterior probability distribution on $K_{\text{sym},0}$

[TABLE]

from which 90% confidence intervals on the curvature of symmetry energy can be extracted. This process is then repeated for $22$ values of chirp mass $\mathcal{M}$ across its feasible range, and then for each interferometer $A$. Appendix C exemplifies this by demonstrating the procedure for one value of chirp mass $\mathcal{M}=1.188\text{ M}_{\odot}$ on interferometer O2, corresponding to GW170817. The results found there are compared to those found in Ref. Carson et al. (2019a) in order to demonstrate the accuracy of our approximated Gaussian $\tilde{\Lambda}$ priors, rather than the full posterior distribution found in Ref. Abbott et al. (2018b). We found that we slightly underestimate the errors in $K_{\text{sym},0}$ by using this method.

There is one important question to analyze here: how do the statistical errors on $K_{\text{sym},0}$ ($\sigma_{A}$ in $P_{A}(\tilde{\Lambda})$ given in Eq. (11) that enters in Eq. (19)) compare to the systematic errors (covariance $\bm{\Sigma}$ in $P(\tilde{\Lambda},K_{\text{sym},0})$ given in Eq. (16) that also enters in Eq. (19))? As more events are observed and the detector sensitivities $S_{n}(f)$ drop, the statistical errors on the measurement of $K_{\text{sym},0}$ approach zero, and the overall errors limit closer to the systematic error “wall” introduced from the EoS variation in the universal relations. We study this effect by first plotting the overall errors on $K_{\text{sym},0}$ as a function of chirp mass, defined to be the one-sided 90% confidence interval on the posterior distribution of $K_{\text{sym},0}$. Following this, we define the systematic errors to be the one-sided 90% confidence interval of $K_{\text{sym},0}$ in the two-dimensional probability distribution evaluated at the central value $\mu_{\tilde{\Lambda}}(\mathcal{M})$ of the $\tilde{\Lambda}$ prior distribution shown in Fig. 6. Equivalently, the fixed diagonal $K_{\text{sym},0}$ coefficient of the Gaussian argument $\exp[-\Sigma^{-1}_{\text{KK}}(K_{\text{sym},0}-\langle K_{\text{sym},0}\rangle)^{2}/2+\dotsc]$ shows the systematic errors to be exactly equal to $(\Sigma_{\text{KK}}^{-1})^{-1/2}$.

Figure 1 displays the results of the above described procedure; plotting the (one-sided 90% confidence interval) overall and systematic errors on the measurement of $K_{\text{sym},0}$ as a function of chirp mass. We observe here the presence of a minimum in the uncertainties with respect to the chirp mass - a relic originating from the correlations between $K_{\text{sym},0}$ and $\tilde{\Lambda}$ seen in Fig. 4, which similarly observe a maximum at the same chirp mass (and thus minimum EoS variation that generates systematic errors). We do note, however, that while previous analyses by Refs. Alam et al. (2016); Malik et al. (2018) required high correlations for the computation of constraints555Refs. Alam et al. (2016); Malik et al. (2018) assumed the relationship between $\Lambda_{1.4}$ (the tidal deformability at $1.4\text{ M}_{\odot}$) and nuclear parameters to lay exactly on the best-fit line between the two. Thus, high degrees of correlation were absolutely necessary for accuracy on this claim., our analysis does not, as all covariances between $\tilde{\Lambda}$ and $K_{\text{sym},0}$ are taken into account by the two-dimensional probability distribution of Eq. (16).

Observe also how, as predicted, the statistical errors drop as the more sensitive detectors are analyzed, reducing to almost zero as the overall errors limit to the fixed systematic error “wall”. The overall errors on the highly-sensitive third generation interferometers CE and ET are indistinguishable from the systematic errors – indicating that the error budget is highly dominated by systematics at this point. Once the errors are dominated by systematics, improving detector sensitivities or observing new events will not aid in the further constraint of $K_{\text{sym},0}$. This indicates the urgent need to reduce the systematic errors found in the EoS-variation of the universal relations for the Voyager-class detectors and beyond666A similar conclusion was reached in Ref. Carson et al. (2019b), where the detector statistical errors became comparable to the systematic errors from the binary Love universal relations for future detectors Voyager and beyond..

IV.2 Multiple events

The future of GW astronomy will become quite busy in terms of detected events. For example, future GW interferometer Cosmic Explorer will be detecting anywhere from $3\times 10^{5}$ to $4\times 10^{6}$ Carson et al. (2019b) binary NS merger events within its horizon distance per year - a staggering number which will certainly help reduce the statistical errors on tidal measurements. How does one account for this effect when studying the uncertainties in future, undetected, events? The dominant tidal parameter in the gravitational waveform, $\tilde{\Lambda}$, depends strongly on the subsequent masses in the binary system, something difficult to predict beforehand. Ultimately, this prevents one from combining the uncertainties on $\tilde{\Lambda}$ for multiple events.

Fortunately, this can be remedied by following in the footsteps of Ref. Carson et al. (2019b), where we reparameterized the gravitational waveform to instead consider the $\Lambda_{1.4}$ and $\Lambda_{1.4}^{\prime}$ tidal coefficients, generated by Taylor expanding the tidal deformability $\Lambda$ about the reference mass of $m_{0}=1.4\text{ M}_{\odot}$ Del Pozzo et al. (2013); Yagi and Yunes (2016):777We note here that a linear truncation of this Taylor series is valid for our purposes. By taking into account an additional quadratic term identified by $\Lambda_{1.4}^{\prime\prime}$, we found a reduction in measurement accuracy in $\Lambda_{1.4}$ by only $\leq 5\%$, across various detectors and values of chirp mass.

[TABLE]

Here, $\Lambda_{1.4}\equiv\Lambda|_{1.4M_{\odot}}$ and $\Lambda_{1.4}^{\prime}\equiv-d\Lambda/d\ln m|_{1.4M_{\odot}}$ are the dimensionless tidal deformability and its slope at $1.4\text{ M}_{\odot}$, and they do not depend on the individual NS masses $m$ or any combination thereof (however they do depend on the fiducial mass value $m_{0}$ chosen). Therefore, they are identical for every future binary NS merger event, and may be combined in uncertainty.

Similar to the correlation between $\tilde{\Lambda}$ and $K_{\mathrm{sym,0}}$, we also find a correlation between $\Lambda_{1.4}$ and $K_{\mathrm{sym,0}}$. We constructed a 2D Gaussian distribution $P(\Lambda_{1.4},K_{\text{sym},0})$ similar to $P(\tilde{\Lambda},K_{\text{sym},0})$ in Eq. (16) and find

[TABLE]

for $\bm{x}=(\Lambda_{1.4},K_{\mathrm{sym,0}})$. Notice that both $\bm{\Sigma}$ and $\bm{\mu}$ are independent of $\mathcal{M}$ in this case.

In this section, we repeat the analysis performed in Sec. IV.1 using the combined uncertainties on $\Lambda_{1.4}$ from $N_{A}$ unique events with chirp mass $1.188\text{ M}_{\odot}$, corresponding to the number of observed binary NS mergers within one observing year on detector $A$. We refer to App. D for details on how to combine information from multiple events, which closely follows App. A of Carson et al. (2019b). Fiducial values of $\Lambda_{1.4}$ and $\Lambda^{\prime}_{1.4}$ were computed to be the mean values of $\Lambda|_{1.4M_{\odot}}$ and $-d\Lambda/d\ln m|_{1.4M_{\odot}}$ from the “LVC constrained” EoSs. Figure 9 shows how the combined-event uncertainties on $K_{\text{sym},0}$ for the fixed chirp mass of $1.188\text{ M}_{\odot}$ further become saturated on the aLIGO, and A+ detectors as well. As was shown in Fig. 1, the single-event uncertainties on $K_{\text{sym},0}$ become dominated by systematic errors for Voyager-class detectors and beyond, and thus there is not much point in stacking multiple events for these detectors to further reduce statistical errors on $\Lambda_{1.4}$.

V Reducing systematic errors via multidimensional correlations

Let us now consider how we can reduce the systematic “walls” present in Fig. 1. In Sec. IV.2, this was computed by evaluating the two-dimensional probability distribution between $K_{\text{sym},0}$ and $\Lambda_{1.4}$ at the fiducial value of $\Lambda_{1.4}$, and then finding the 90% confidence interval of the resulting probability distribution of $K_{\text{sym},0}$ to yield $\sim 104$ MeV. We here construct multidimensional correlations among $K_{\mathrm{sym,0}}$ and $\Lambda_{m_{x}}$ at a few different masses $m_{x}$ (since we expect to detect GWs from binary NSs with different masses with future observations) to see how adding information of the tidal deformability at multiple different masses may help us to reduce the systematic errors on $K_{\mathrm{sym,0}}$.

Let us begin by using $\Lambda$ at two different masses $m_{x}$ and $m_{y}$. This requires us to find a three-dimensional correlation among $K_{\text{sym},0}$, $\Lambda_{x}(\equiv\Lambda_{m_{x}})$, and $\Lambda_{y}(\equiv\Lambda_{m_{y}})$, and construct a three-dimensional Gaussian distribution $P(K_{\text{sym},0},\Lambda_{x},\Lambda_{y})$. Figure 10 shows an example of such a distribution for the case of $m_{x}=1.3\text{ M}_{\odot}$ and $m_{y}=1.6\text{ M}_{\odot}$. The systematic error is then computed by evaluating the three-dimensional distribution at the fiducial values of $\Lambda_{x}$ and $\Lambda_{y}$, and then evaluating the resulting one-dimensional $K_{\text{sym},0}$ probability distribution at the 90% confidence interval.

Figure 11 displays the resulting systematic uncertainties on $K_{\text{sym},0}$ using canonical masses $m_{x}$ and $m_{y}$ between $1\text{ M}_{\odot}$ and $2\text{ M}_{\odot}$. Observe that the systematic errors can be reduced by setting both $m_{x}$ and $m_{y}$ to be large or small. However, this means that $m_{x}\approx m_{y}$, which corresponds to effectively using a two dimensional probability distribution. Thus, in this case, having additional information on $\Lambda$ at a different mass does not help to reduce the systematic errors.

It may sound strange that adding more pieces of information does not help to reduce the systematic errors. Let us explain why this is the case by comparing the systematic errors at $(m_{x},m_{y})=(1,2)M_{\odot}$ and $(m_{x},m_{y})=(2,2)M_{\odot}$. Figure 12 compares the two-dimensional 90% contours between $K_{\text{sym},0}$ and $\Lambda_{2.0}$ from two different methods by computing (i) directly the two-dimensional probability distribution from Eq. (16), and (ii) the three-dimensional probability distribution between $K_{\text{sym},0}$, $\Lambda_{2.0}$, and $\Lambda_{1.0}$, and then evaluating it at the fiducial value of $\Lambda_{1.0}$. We observe that while the contour from the first case has a larger area (and value of $|\bm{\Sigma}|$) as expected due to the use of less information, it becomes distorted such that the systematic uncertainty (along the dashed horizontal line corresponding to the fiducial value of $\Lambda_{2.0}$) becomes smaller than that from the first case.

Let us now consider using $\Lambda$ at three different masses, $m_{x}$, $m_{y}$ and $m_{z}$. This requires us to find a four-dimensional correlation and construct the four-dimensional Gaussian probability distribution $P(K_{\text{sym},0},\Lambda_{x},\Lambda_{y},\Lambda_{z})$. We fix $m_{z}=1.5\text{ M}_{\odot}$, and allow $m_{x}$ and $m_{y}$ to vary between $[1.0,2.0]\text{ M}_{\odot}$. Similar to the process used previously, this probability distribution is evaluated at the fiducial values of $\Lambda_{x}$, $\Lambda_{y}$ and $\Lambda_{1.5}$:

[TABLE]

The resulting 90% confidence intervals are presented in Fig. 2 for the entire range of $m_{x}$ and $m_{y}$ mass values. We observe that by including information about binaries with large, medium, and small masses together, the systematic errors can be improved drastically, down to $\sim 74$ MeV. We also see that along the diagonal line of $m_{x}=m_{y}$, the four-dimensional probability distribution $P(K_{\text{sym},0},\Lambda_{x},\Lambda_{x},\Lambda_{1.5})$ reduces to the three-dimensional case $P(K_{\text{sym},0},\Lambda_{x},\Lambda_{1.5})$, with uncertainties $\sim 104$ MeV approaching that of Fig. 11 along the horizontal dashed line.

VI Discussions

Are there any other ways to further improve the constraints on $K_{\text{sym},0}$ using observed GW events? One might think that the constraint on $\tilde{\Lambda}$ with GW170817 may help in this direction. However, the restriction of data in only the $\tilde{\Lambda}$ dimension does not help as systematic errors are found by evaluating the scattering width in the $K_{\text{sym},0}$ direction.

Finally, we briefly discuss the possibility of NSs with strong first-order phase transitions from hadronic to quark matter in the core, as described in Ref. Paschalidis et al. (2018). With high enough observed chirp masses $\mathcal{M}$, future binary NS merger events could potentially be composed of one or both hybrid stars (HSs) with quark-matter cores. The tidal deformabilities and thus, the nuclear parameters, depend on such structure, and could potentially disagree between events with varying chirp masses and combinations of NS/HS Montana et al. (2019). Thus, significant variations between nuclear parameter measurements with future GW observations with varying chirp masses could potentially present evidence of strong phase transitions at around $2-3$ times the nuclear saturation density. If such transitions are present at sufficiently low-densities, then the nuclear matter parameters will be further decoupled from the GW observations and thus our lower limit for the nuclear matter parameter uncertainties will increase. Alternatively, similar measurements of nuclear parameters could either indicate a pure hadronic matter EoS, or phase transitions occurring at higher nuclear densities. The structure of such high-density transitions could be probed by the GW post-merger oscillation signal.

As we showed in the previous section, one can use multidimensional correlations to reduce the systematic errors. Instead of using tidal deformabilities from different NS masses obtained from GWs alone, one can consider combining information from multi-messenger observations. For example, pulse profiling techniques by NASA’s Neutron star Interior Composition Explorer (NICER) instrument may provide high-precision measurements on the NS radius down to $5\%$ Ozel et al. (2016); Sieniawska et al. (2018). Thus, one can construct multidimensional correlations among nuclear parameters, tidal deformabilities and radii of NSs. The work along this direction is currently in progress Schumacher et al. (tion).

Acknowledgments

We thank David Nichols for his illuminating advice on conditional probability distributions. Z.C. and K.Y. acknowledge 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 Original versus restricted sets of EoS

In this appendix, we show how the restriction of EoSs described in Sec. II.3 used in the current analysis impacts our observations, as compared to the original set of EoSs used in Ref. Carson et al. (2019a). Figure 13 shows a comparison between the two-dimensional probability distributions $P(K_{\text{sym},0},\tilde{\Lambda})$ resulting from each set of EoSs. We see that while the restriction to EoSs does indeed shrink the 90% confidence intervals in the direction of correlation, the widths are approximately equal at the fiducial value of $\tilde{\Lambda}$ (where the systematic errors are analyzed). This shows that while using a subset of EoSs may reduce the overall error888We indeed observed large reductions in the overall errors found in Fig. 1 when using the restricted EoSs rather than the original ones, while the level of systematics stayed constant at $\sim 104$ MeV., the level of systematic errors will remain mostly fixed. Our result is consistent, for example, with Ref. Carriere et al. (2003), which found that the correlation between $L_{0}$ and the radius of a $1.4\text{ M}_{\odot}$ neutron star was weak because of the contribution of the high-density component of the EoS.

Appendix B Multiplicative combinations of nuclear parameters

In this appendix, we discuss the feasibility of using multiplicative combinations of nuclear parameters, such as $K_{0}L_{0}^{\eta}$, rather than the linear combinations such as $K_{0}+\alpha L_{0}$ considered in Refs. Alam et al. (2016); Malik et al. (2018); Carson et al. (2019a). Here we consider the following multiplicative combinations for comparison purposes: $K_{0}L_{0}^{\eta}$, $M_{0}L_{0}^{\mu}$, and $K_{\text{sym},0}L_{0}^{\nu}$, where coefficients $\eta$, $\mu$, and $\nu$ are similarly chosen to achieve maximal correlation. Such multiplicative combinations are similar to those considered in Sotani et al. (2014); Silva et al. (2016).

Figure 14 presents the correlations between $\tilde{\Lambda}$ and all 6 multiplicative and linear combinations of nuclear parameters considered in this analysis. While the two classes of nuclear parameter combinations produce very similar correlations with $\tilde{\Lambda}$, we observe that the linear cases slightly outperform the multiplicative cases for nearly all values of chirp mass. Similarly, repeating the analysis999Because the two-dimensional probability distribution is now between $\tilde{\Lambda}$ and $K_{\text{sym},0}L_{0}^{\nu}$, an additional marginalization over $L_{0}$: $\int^{\infty}_{-\infty}P(K_{\text{sym},0}L_{0}^{\nu})P(L_{0})dL_{0}$ must be performed ($P(L_{0})$ is an additional prior distribution on $L_{0}$ given by Refs. Lattimer and Lim (2013); Lattimer and Steiner (2014); Tews et al. (2017); Oertel et al. (2017)) in order to extract the posterior distribution on $K_{\text{sym},0}$. found in Sec. IV returns constraints on $K_{\text{sym},0}$ to be slightly worse than that considered in the the main analysis, due to the additional inclusion of uncertainties from nuclear parameter $L_{0}$. One arrives at a similar conclusion if one uses a linear combination with $\gamma\neq 0$ Carson et al. (2019a).

From this evidence, we conclude with the remarks that the multiplicative combinations of nuclear parameters offer nothing new in terms of enhanced constraints on nuclear parameters. The multiplicative combinations of nuclear parameters slightly under-perform their linear combination counterparts in terms of correlations with $\tilde{\Lambda}$. Thus, we neglect their use and continue our analysis as was done previously in Carson et al. (2019a).

Appendix C Example computation of the $K_{\text{sym},0}$ posterior distribution

In this appendix, we demonstrate the process of computing the posterior distribution on $K_{\text{sym},0}$ (used in Sec. IV) for one value of chirp mass, $\mathcal{M}=1.188\text{ M}_{\odot}$, corresponding to GW170817 on detector O2. This case corresponds to the large dot in Fig. 1. Referring to Figs. 6 and 7, we observe that the mean and root-mean-square $\tilde{\Lambda}$ values for O2 detector sensitivity at $\mathcal{M}=1.188\text{ M}_{\odot}$ are given by $\mu_{\tilde{\Lambda}}=430.8$ and $\sigma_{\text{O2}}=172.5$, respectively. This results in a prior distribution on $\tilde{\Lambda}$ shown in Fig. 15, given by

[TABLE]

We additionally show the true posterior distribution on $\tilde{\Lambda}$ derived in Ref. Abbott et al. (2018b), which was used as a prior in our original analysis found in Ref. Carson et al. (2019a) for comparison purposes.

Following along with Sec. IV.1, we generate the two-dimensional probability distribution between $K_{\text{sym},0}$ and $\tilde{\Lambda}$, given by Eq. (16). We find the covariance matrix and the mean vector to be

[TABLE]

for $\bm{x}=(\tilde{\Lambda},K_{\mathrm{sym,0}})$ This results in the two-dimensional probability distribution between $K_{\text{sym},0}$ and $\tilde{\Lambda}$ shown in Fig. 16. The systematic error on $K_{\text{sym},0}$ is then computed by evaluating the 90% confidence interval width of the distribution in the $K_{\text{sym},0}$-dimension at $\mu_{\tilde{\Lambda}}=430.8$, corresponding to the mean of the prior distribution in $\tilde{\Lambda}$. The resulting (one-sided 90% confidence level) systematic errors for this case are found to be $\sigma_{\text{sys}}=104.6\text{ MeV}$.

Next, we obtain the one-dimensional conditional probability distributions on $K_{\text{sym,0}}$ given tidal deformability observations of $\tilde{\Lambda}_{\text{obs}}$. By following Eq. (18), this is simply given by

[TABLE]

Finally, the posterior distribution on $K_{\text{sym},0}$ can be computed by combining the one-dimensional conditional probability distributions with the prior distribution on $\tilde{\Lambda}$, and then integrating over all observations of $\tilde{\Lambda}$:

[TABLE]

Figure 17 finally displays the resulting posterior distribution on $K_{\text{sym},0}$, with a mean value of $-115^{+75}_{-73}\text{ MeV}$, giving a 90% confidence interval of $-227\text{ MeV}\leq K_{\text{sym},0}\leq 7\text{ MeV}$ (or a one-sided $90\%$ confidence interval of $117$ MeV corresponding to the maroon dot in Fig. 1). Comparing this to the resulting posterior distribution (also shown in Fig. 17) found in Ref. Carson et al. (2019a) giving a 90% confidence interval of $-285\text{ MeV}\leq K_{\text{sym},0}\leq 7\text{ MeV}$, we find that this approximation of $\tilde{\Lambda}$ priors slightly underestimates the errors in $K_{\text{sym},0}$, but otherwise works quite well. We also see that this approximation of $\tilde{\Lambda}$ priors skews the distribution less so than the LVC posterior, making it more normally distributed. We also note that here, we utilize a restricted set of EoSs compared to that done in Ref. Carson et al. (2019a), resulting in a slightly more accurate posterior distribution on $K_{\text{sym},0}$.

Appendix D Combining uncertainties from multiple events

In this appendix, we explain how one can combine statistical uncertainties on certain parameters from multiple events. We closely follow App. A of Carson et al. (2019b). To estimate the number of detected binary NS merger events $N_{A}$ on interferometer $A$, we integrate the local binary NS merger rate across all redshift values within detector $A$’s horizon redshift $z_{h}$ with SNR thresholds of $\rho_{\text{th}}=8$. Following the process used in Ref. Carson et al. (2019b), this is given by

[TABLE]

In our cosmology, we choose $a_{0}r_{1}(z)$, $\frac{d\tau}{dz}$, and $r(z)$ to be

[TABLE]

[TABLE]

[TABLE]

where $H_{0}=70\text{ km s}^{-1}\text{ Mpc}^{-1}$ is the local Hubble constant, $\Omega_{\Lambda}=0.67$ is the universe’s vacuum energy density, $\mathcal{R}=1540^{+3200}_{-1220}\text{ Gpc}^{-3}\text{ yr}^{-1}$ is the local binary NS coalescence rate density Abbott et al. (2017b), and $\Delta\tau_{0}$ is chosen to be a one year observation period. Doing so gives the rates found in Table VIII of Ref. Carson et al. (2019b). For demonstration purposes, we choose to use the lower and upper limits of these rates which are found to be $(2.0\times 10^{0},3.0\times 10^{2})$, $(1.6\times 10^{2},2.4\times 10^{3})$, $(2.2\times 10^{3},3.2\times 10^{4})$, $(7.2\times 10^{4},1.1\times 10^{6})$, and $(3.0\times 10^{5},4.4\times 10^{6})$ events per year for aLIGO, A+, Voyager, CE, and ET respectively.

Similarly following Ref. Carson et al. (2019b), the combined uncertainty $\sigma_{N_{A}}$ is computed by first simulating a population of $N_{A}$ events from the SNR ($\rho$) probability distribution: Schutz (2011); Chen and Holz (2014)

[TABLE]

with an SNR threshold of $\rho_{\text{th}}=8$. The combined population root-mean-square error on $\Lambda_{1.4}$ is then computed by integrating over all $N_{A}$ sources at various redshifts $z$:

[TABLE]

Here, $\sigma^{A}(z)$ is the redshift dependence of the root-mean-square error on detector $A$, evaluated via Fisher analyses at various redshifts up to $z_{h}$.