Constraining the neutronmatter equation of state with gravitational waves
Abstract
We show how observations of \glsGWs from \glsBNS mergers over the next few years can be combined with insights from nuclear physics to obtain useful constraints on the \glsEoS of dense matter, in particular, constraining the neutronmatter \glsEoS to within 20% between one and two times the nuclear saturation density . Using Fisher information methods, we combine observational constraints from simulated \glsBNS merger events drawn from various population models with independent measurements of the neutron star radii expected from xray astronomy (the \glsNICER observations in particular) to directly constrain nuclear physics parameters. To parameterize the nuclear \glsEoS, we use a different approach, expanding from pure nuclear matter rather than from symmetric nuclear matter to make use of recent \glsQMC calculations. This method eschews the need to invoke the socalled parabolic approximation to extrapolate from symmetric nuclear matter, allowing us to directly constrain the neutronmatter \glsEoS. Using a principal component analysis, we identify the combination of parameters most tightly constrained by observational data. We discuss sensitivity to various effects such as different component masses through populationmodel sensitivity, phase transitions in the core \glsEoS, and large deviations from the central parameter values.
pacs:
Contents
I Introduction
The detection of gravitational waves from the binary neutron star (\glsBNS) merger \glsGW170817 by the \glsaLIGO detectors Aasi et al. (2015) in Hanford, \myscwa (\mysclho) and Livingston, \myscla (\myscllo) and the \glsVirgo detector Acernese et al. (2015) ushered in the era of multimessenger astronomy with gravitational waves Abbott et al. (2017a, b). This has been instrumental in launching novel ways of constraining cosmological parameters Abbott et al. (2017c); Chen et al. (2018); Nair et al. (2018); SoaresSantos et al. (2019), on the one hand, and neutron star equation of state (\glsEoS) parameters, on the other hand Abbott et al. (2017a, 2018a). In a \glsBNS system the \glsNS masses and their \glsEoS determine how much quadrupolar deformation their tidal fields are able to induce in each other. The two are related by the tidal deformability parameter as . It is now well understood that the tidal deformability parameters of both \glsplNS in a double neutron star system affect the phase of the gravitational wave signal during the late stages of the inspiral Flanagan and Hinderer (2008).
Recent articles that followed discovery of \glsGW170817 have shown that upper bounds on the dimensionless tidal deformability of the neutron stars obtained from gravitational wave data analysis provide constraints on the \glsEoS of dense matter encountered inside neutron stars De et al. (2018); Tews et al. (2018a); Abbott et al. (2019). This is a great opportunity and challenge for several reasons: neutron rich matter, although relevant for many applications, is not easily accessible in experiments, while theoretical approaches require solving the difficult quantum manybody problem and lack a precise characterization of the underlying interactions. Observational constraints provide an anchor for nuclear theory in this uncertain regime, allowing one to extrapolate lowdensity and symmetric properties of nuclear matter to significantly improve constraints on neutronrich matter at higher densities.
In this article we discuss how we can extract more detailed information about the properties of dense neutronrich matter and neutron stars during the next few years with more \glsGW detections and measurements of neutron star radii expected from xray astronomy, and highlight the importance of an informed parameterization of the dense matter \glsEoS. We make the reasonable assumption that all neutron stars are described by the same \glsEoS. Further, modern nuclear Hamiltonians based on chiral effective field theory provide a systematic momentum expansion of two and manybody nuclear forces. This, combined with advanced computational methods to solve the nonrelativistic quantum manybody problem, now allows us to calculate the \glsEoS of pure neutron matter up to nucleon number density , where nucleons per is the average nucleon density inside large nuclei (corresponding to a mass density ) Tews et al. (2018b). Interestingly, there is a convergence of different ab initio methods based on realistic microscopic Hamiltonians that account for two and three neutron forces Gandolfi et al. (2009); *Gandolfi:2010b; *Gandolfi:2012; *Gandolfi:2014a. These calculations suggest that the functional form of the \glsEoS of pure neutron in the density interval to is well determined. We use this information to parameterize the \glsEoS and show how it helps with the analysis of multiple \glsBNS detections and provide tighter and more useful constraints for dense matter physics. In turn, these constraints for the \glsEoS of pure neutron in the density interval where calculations are feasible will provide new insights for nuclear physics.
Our study differs from earlier work in the following aspects:

We incorporate insights about neutronrich matter obtained from nuclear physics by implementing a new parameterization of nuclear equation of state and identify parameters that can be best constrained by \glsGW observations.

We quantify how constraints on these parameters and on the pressure of neutron matter in the density interval to will improve with the number of detections.

Our analysis uses a numerical relativity based tidal waveform model.

We study the effect of different population synthesis models on the accuracy with which \glsEoS parameters can be measured with \glsplGW and use several thousand binary neutron star source simulations to assess errors in \glsEoS parameter measurements.

While a nearby event like \glsGW170817 at \glsaLIGO design sensitivity would significantly constrain the properties of neutron matter, we show that similar constraints can be obtained from about 15 events beyond .
We begin with a summary of our results in section II, then describe how we have parameterized the dense matter \glsEoS in section III. In section IV we discuss how we obtain constraints from the gravitational waveform of simulated merger events. Finally, we discuss details of the method we use to obtain these constrains in section VI.
Ii Results
Our main result is that even a handful of \glsGW observations of \glsBNS mergers will provide the most stringent constraints on the lowtemperature equation of state of dense neutron matter in the density interval between . This is summarized in fig. 1, which shows how the constraints on the pressure of pure neutron matter improve as a function of additional \glsNICER or \glsLIGO observations. We start from the errors listed in table 1, which, for the purposes of this analysis, we interpret as uncorrelated normal errors for the parameters. This gives the upper dotted line labeled “Nuclear”.
To this, we add the following constraints:

Constraints from a simulated binary with similar masses and distance to \glsGW170817 but at \glsaLIGO design sensitivity.

\Gls
GW observations at \glsaLIGO design sensitivity of distant simulated merger events from population model \glsCompleteStdAsubsolarNSNS as described in section IV. To estimate the variance possible within the population model, we sample 500 different populations, each containing , and plot the (68th percentile) error bands as shaded regions.
This analysis demonstrates several key points: A nearby event such as \glsGW170817 is comparable to a dozen or so events from . The \glsNICER constraints are comparable to a single \glsLIGO observation from a distant population sample having low \glsSNR, however, nearby or multiple accumulated \glsLIGO events yield significant improvement. After about observation, we observe rather limited improvement from additional . This can also be seen in fig. 2, which shows how the constraints improve as a function of the number of observations.
One caveat: these constraints assume Gaussian errors and linear error propagation. A proper analysis requires a much more expensive Bayesian approach (see, e.g., Ref. Agathos et al. (2015)). To assess the nonlinear effects, we provide similar plots for comparison in EPA () for the different central values listed in table 2.
To put these results in perspective, consider the nuclear symmetry energy and the slope of its density dependence , \cref@addtoresetequationparentequation
(1a)  
(1b) 
where is the energyperparticle of uniform nuclear matter. If the socalled parabolic approximation holds at saturation ( – see eq. 5 and the surrounding discussion), then upcoming neutron skin experiments Horowitz et al. (2012); *Horowitz:2014; *Horowitz:2014a expect to constrain with a possible reduction to with a followup experiment. This is comparable to combined constraints from ab initio calculations Hebeler and Schwenk (2010); *Hebeler:2013; Wlazłowski et al. (2014); Gandolfi et al. (2014b); Lynn et al. (2015) and astrophysical observations Page and Reddy (2006); Gandolfi et al. (2009); *Gandolfi:2010b; *Gandolfi:2012; *Gandolfi:2014a; Steiner and Gandolfi (2012); *Steiner:2013; Lattimer and Steiner (2014). From our analysis we thus see that \glsGW observations alone could have an impact at the level corresponding to .
Iii Parameterization of the Nuclear Equation of State
\glsresetEoS To relate the nuclear equation of state to the structure of neutron stars, we must first characterized the \glsEoS of nuclear matter. This is conveniently parameterized by the energy density as a function of the baryon number density , which is the sum of the neutron and proton number densities. Simple approximation for this function in terms of polytropes are often a starting point for astrophysical analysis. Indeed, many families of nuclear \glsEoS can be characterized quite well by a simple set of piecewise polytropes Read et al. (2009).
Our approach here, however, is to directly express in terms of nuclear physics parameters. This approach allows one to directly assess how observations translate into constraints on nuclear physics. We shall demonstrate this by providing constraints on the pressure of pure neutron matter , which is inaccessible from a general polytropic analysis (fig. 1).
It is useful to divide the neutron star interior into four regions: the outer crust, the inner crust, the outer core, and the inner core. The radial extent of the outer crust, which is composed neutronrich nuclei embedded in a electron gas, is only a few hundred meters and its contribution to the neutron star mass is negligible. The \glsEoS of the outer crust is well understood and depends weakly on the composition of nuclei present. The inner crust extends from to , has radial thickness , and contains a modest fraction of the mass. Here, exotic neutronrich nuclei are embedded in a dense liquid of neutrons and electrons, as described by the \glsCLDM in section III.1. The outer core is a liquid composed primarily of neutrons and a small (few percent) admixture of protons, electrons, and muons. It extends from to where the description of matter in terms of nucleons interacting with static potentials is expected to break down. The inner core extends to higher densities, and we switch here to the speedofsound parameterization discussed in section III.3.
On dimensional grounds one expects the dimensionless tidal deformability to be related to with for Nelson et al. (2018). Although the \glsEoS around intermediate densities dominates the 7th moment of energy distribution for massive neutron stars, the inner crust also makes a large contribution to for lowmass stars (which are believed to be more common in binary neutron star systems). This contribution is shown by the dashed (purple) lines at the bottom of fig. 3. Thus, it is important to provide a unified description of the \glsEoS of the inner crust and the outer core in any analysis that aims to constrain the \glsEoS using \glsGW observations of binary neutron stars.
iii.1 Compressible Liquid Drop Model
\glsresetCLDM The \glsCLDM (see Haensel et al. (2007); *Chamel:2008) provides a unified \glsEoS connecting a fixed outer crust for (for which we use the data in Table 4 of Sharma et al. (2015)) to the inner core \glsEoS. In the inner crust, the \glsCLDM constructs spherical nuclei in a spherical WignerSeitz cell, ensuring equilibrium with surrounding neutron and lepton gases by establishing both electric and equilibrium. This is similar to the approach taken in Fortin et al. (2016); Zdunik et al. (2016), but differs in how we define the nuclear matter \glsEoS . Instead of using obtained from specific models based on effective Hamiltonians solved in the mean field approximation to reproduce empirical parameters like nuclear saturation properties, we use what we believe is close to a minimal phenomenological parameterization that directly encodes properties that can either be measured or calculated reliably. The advantage of our approach is that these parameters are directly connected with the unified \glsEoS, allowing us to provide a full covariance analysis linking nuclear parameters with neutron star observables.
Although the use of a spherical WignerSeitz cell precludes the possibility of pasta phases Ravenhall et al. (1983) the errors incurred by the WignerSeitz approximation for different lattice structures are less then 0.5% (see e.g. Chamel et al. (2007); Chamel and Haensel (2008)).
Our implementation of the \glsCLDM introduces two effective parameters: the surface tension and the parameter which characterizes the isospin dependence of the surface tension Lattimer et al. (1985) (see EPA () for the exact form used), where is the isospin asymmetry. We fix the parameter to smoothly match the tabulated outer crust equation of state, leaving free the single parameter . Additionally, we include as a parameter a suppression factor for the Coulomb interaction to allow for the diffusivity of the proton charge distribution (see the discussion in Steiner (2012)). As will be shown in section II, these parameters have negligible effects on the constructed \glsEOS.
This approach allows for a small firstorder phase transition from the region modeled by the \glsCLDM to homogeneous nuclear matter. With our parameters, this phase transition is weak: .
To establish equilibrium we include leptons modeled as a Fermi gas of electrons (and muons at sufficiently high densities) in the \glsTF approximation.
\GlsentrynameCLDM parameters:

Neutron Matter  Inner Core  
\glsEoS  []  []  []  
\glsCentral2c  
\glsSoft2c  
\glsStiff2c  
\glsSoft2c_Stiff  
\glsStiff2c_Soft  
\glsCentral2c_low_Ec  
\glsCentral2c_high_Ec  
\glsCentral2c_low_C_max  
[]  
\glsCentral2c_trans 
iii.2 Homogeneous Nuclear Matter
One of the main new features of our analysis is to parameterize the nuclearmatter \glsEoS as an expansion in the proton fraction from pure neutron matter to symmetric neutron matter. This is in contrast to the common approach of expanding about symmetric nuclear matter in powers of the isospin asymmetry . The common approach allows one to directly connect experimentally relevant properties of symmetric nuclear matter to properties of neutron matter. This connection, however, is generally predicated on the socalled parabolic approximation, which is valid only if quadratic terms dominate over quartic and higherorder terms. While there is some support for this below saturation density from relativistic \glsDBHF calculations Lee et al. (1998), Gogny forces GonzalezBoquera et al. (2017), and other perturbative techniques (see Li et al. (2008) for a review), it is not well established at higher densities. Indeed virtually any form of neutronmatter \glsEoS can be accommodated with quartic terms without spoiling global mass fits Bulgac et al. (2018). For this reason, we start with a parameterization of pure neutron matter, then use the properties of symmetric nuclear matter to constrain the extrapolation in the proton fraction .
To describe pure neutron matter we use a double polytrope for the energy per particle:
(2) 
where is the nucleon mass, is a constant (approximately the nuclear saturation density), and , , , and are four \glsEoS parameters. This form was found to accurately fit \glsQMC calculations of the \glsEoS using nuclear Hamiltonians with realistic two and threebody forces Gandolfi et al. (2009); *Gandolfi:2010b; *Gandolfi:2012; *Gandolfi:2014a, and is consistent with recent \glsQMC results based on chiral \glsEFT interactions Wlazłowski et al. (2014); Gandolfi et al. (2014b, 2015); Lynn et al. (2015). For small proton fractions , we perform an expansion:
(3) 
where is the proton effective mass, and describes the selfenergy of the proton polaron. This function is presently poorly constrained by \glsQMC and experimental data and all known results are consistent with a simple twoparameter quadratic expansion:
(4) 
where and where returns to zero. (We expect to curve up for higher densities due to the repulsive nature of nuclear threebody interactions).
The additional powers are chosen to match the properties of nuclear matter to quadratic order in the isospin asymmetry and expansion away from saturation :
(5) 
Fitting two even powers, and , and the lack of odd powers uniquely defines the functions through , completing our characterization of the nuclear equation of state in terms of the nuclear saturation density , energy , and incompressibility ; the symmetry energy , slope and incompressibility . Note that a term proportional to is allowed in eq. 5, but our \glsEoS is unconstrained by this term, i.e., does not rely on the parabolic approximation eq. 5.
iii.3 Speed of Sound Parameterization of the Inner Core
Above densities the \glsEoS is virtually unconstrained. The typical approximation at high density is in terms of a polytrope, but we choose a more physically motivated highdensity \glsEoS parameterized in terms of the square of the speed of sound: which approaches the \glsPQCD result at asymptotic densities. Although the form of the function is unknown at finite density, its qualitative form at finite temperature suggests that it may first peak before returning to the asymptotic value Alford (); Tews et al. (2018b). We thus include a simple parameterization as a quadratic polynomial smoothly connecting to the homogeneous equation of state at a fixed transition energy density reaching a maximum at an energy density , then returning to at which it remains for higher densities. This core \glsEoS thus introduces three parameters , , and . To better understand the sensitivity of our results to the properties of the core, we include one slightly different form \glsCentral2c_trans which has a firstorder phase transition with discontinuity at .
iii.4 Parameters
Our equation of state is thus characterized by 18 parameters: and , (\glsCLDM), , , , (symmetric nuclear matter), , , , (symmetry energy), , , , , (neutron matter) , , , (proton polaron), and , , (core). We explore various ranges of these parameters centered about the values listed in table 1, which defines our base \glsCentral2c \glsEoS model. In addition to these central values, we repeat our analysis at a handful of different parameter values, defining the models listed in table 2. Some of these are referred to in the text, but a complete comparison is present in the supplement EPA (). We now discuss how these constraints are derived from gravitational wave observations.
Iv Gravitational Waveform
Gravitational waves from merging binary neutron star systems carry information about the nuclear equations of state. During late stages of inspiral tidal interactions between neutron stars can leave imprints on the \glsGW signal that is otherwise dominated by pointmass contributions. As mentioned earlier, tidal responses of neutron stars can be quantified by the dimensionless tidal deformability parameter , where the second Love number is weakly sensitive to the matter distribution inside the star Flanagan and Hinderer (2008). The strong dependence of on the radius of neutron star allows us to extract information regarding nuclear \glsEoS. Indeed, \glspN theory is able to quantitatively describe the effect of the \glsNS \glsEoS on the signal by parameterizing the waveform in terms of and of component stars Flanagan and Hinderer (2008); Vines et al. (2011).
Gravitational wave observations of inspiraling compact binaries involving neutron stars can therefore constrain Abbott et al. (2017a, 2018a). However, since the constraint on from a single \glsBNS is weak for small to medium \glsSNR events, multiple observations of such systems will be required for remote sources to reduce the statistical error in s and s in order to discern the effects of similar \glsEoS Del Pozzo et al. (2013); Agathos et al. (2015); Bose et al. (2018). Fortunately, tenstohundreds of binaries of this type Abbott et al. (2018b) are expected to be observed over the next several years by the advanced (or “second generation”) \glsLIGO.
We consider only nonspinning neutron stars here because astrophysically their spins are expected to be small when in a \glsBNS system; in particular it is believed that the dimensionless spin parameter Stovall et al. (2018); Abbott et al. (2017a) We plan to study the effect of spin in a future follow up study.
The \glsGW signal from a \glsBNS system in a detector can be expressed as the strain
(6) 
where and denote its amplitude and phase in the time domain. For \glsFIMbased parameter estimation, we work with the Fourier transform of the strain above. This is constructed by adding to the pointparticle part of the TaylorF2 model at 3.5\glspN Buonanno et al. (2009), a phase correction that is taken here to be the Fourier domain tidal waveform, with Padé fits, as prescribed in Dietrich et al. Dietrich et al. (2017).
V Population Models
We employ different sets of stellar evolution model parameters of \glsZAMS binary stars each of which would lead to a binary neutron star system that merges within Hubble time. The differences among stellar evolution models can be large, resulting in appreciable variation in the component mass distribution. Since the tidal deformability parameter is sensitive to the masses, we explore four cases of mass distributions produced by population synthesis studies Dominik et al. (2012). These are more realistic than the uniform or Gaussian distributions owing to the application of stellar evolution mechanism of binary stars including two important factors, namely, metallicity and the nature of the common envelop interaction in the binary.
Metallicity plays the most dominant role in determining the strength of stellar winds in main sequence stars. The larger the metallicity the larger the stellar winds, due to increased scattering crosssection of the electrons. This results in increased mass loss; therefore, the remnant mass left behind at the end of main sequence phase is reduced. This decreases the total baryonic mass content of the supernova engine at the onset of the explosion. In our study, we consider two different variants of metallicities produced by Dominik et al. (2012). In the first case, the stellar evolution model was used with metallicity abundances being the same as solar metallicity, while in the second case 1/10th of solar metallicity was used. The latter is termed to be of subsolar metallicity. Component masses are narrowly peaked for solar metallicity systems while subsolar metalicity system produce a wider mass distribution.
The second most important effect that can change the component masses of \glsBNS systems is the way mass transfer takes place during the common envelop phase of stellar evolution of the binary stars. The mass transfer in the common envelop stage depends on the evolutionary phase of the two stars. In one extreme case, for example, if the common envelop phase is initiated by the star in the Hertzsprung gap stage, it is likely to transfer a significant amount of orbital angular momentum to the entire binary system. This case is denoted by “submodel A” in Dominik et al. (2012). On the other hand, depending on the nature of interaction between the core and the envelop, one possible outcome is that during each common envelop stage for Hertzsprung gap donor stars the outer envelope acquires the significant part of the orbital angular momentum and gets ejected from the system, leaving behind the cores of the two stars to inspiral. This case is denoted by “submodel B” in Dominik et al. (2012). Furthermore, a higher metallicity in the parent star can result in greater mass loss and consequently a less massive remnant. Therefore, we employ \glsNS populations resulting from solar metallicity stars as well as those with 10% of solar metallicity. These different characteristics lead to the following four categories of population models studied here: \glsresetStdASolar \glsresetCompleteStdAsubsolarNSNS \glsresetCompleteStdBsolarNSNS \glsresetCompleteStdBsubsolarNSNS \glsresetUniform
 \GlsStdASolar

These are binary \glsNS populations produced by solar metallicity stars of the submodel A type.
 \GlsCompleteStdAsubsolarNSNS

These are binary \glsNS populations produced by subsolar metallicity stars of the submodel A type.
 \GlsCompleteStdBsolarNSNS

These are binary \glsNS populations produced by solar metallicity stars of the submodel B type.
 \GlsCompleteStdBsubsolarNSNS

These are binary \glsNS populations produced by subsolar metallicity stars of the submodel B type.
 \GlsUniform

Uniform sampling of neutron stars with masses between and .
Vi Statistics and Methods
Given a particular parameterization of the \glsEoS, we compute the mass , radius , and tidal deformability parameter of a \glsNS with a given central density by solving the \glsTOV equations (see e.g. Hinderer (2008); Postnikov et al. (2010)). The signals (gravitational waveforms) from merging neutron stars is computed with the numerical relativity based frequencydomain model Dietrich et al. (2017) mentioned above. From those waveforms, we compute the corresponding \glsFIM characterizing the correlated uncertainties of the masses, and , and the tidal deformabilities, and (maximizing the matchedfilter over the source distance, signal time, and phase at colaescence Ajith and Bose (2009)), to estimate the information obtainable in a merger event at \glsaLIGO design sensitivity, as described below.
Statistical Analysis
To estimate how large the noiselimited errors are of the \glsBNS parameters , we begin by modeling the measured values after the \glsplMLE Helstrom (1995). Owing to noise, the \glsMLE will fluctuate about the respective true values, i.e., , where is the random error. The extent of these fluctuations is estimated by the elements of the variancecovariance matrix, Helstrom (1995).
The matrix is bounded by the signal via the CramerRao inequality, which states that
(7) 
where is the \glsFIM:
(8) 
Above, is the partial derivative with respect to the parameter and is the onesided noise \glsPSD Helstrom (1995). We take the latter to be the zerodetuned highpower \glsPSD for \glsaLIGO aLI (2010). Therefore, gives the lower bound on the \glsrms error in estimating . The two are equal in the limit of large \glsSNR (see, e.g., Vallisneri (2008)). The error estimates listed here are the obtained from the \glsFIM.
The \glsFIM method is known to underestimate the error in the estimation of the masses Rodriguez et al. (2013). We therefore used errorestimates for totalmass and massratio (i.e., the ratio of the lighter mass to heavier mass) that were obtained with Bayesian methods in Ref. Rodriguez et al. (2014), and set them such that the error is and , respectively, at a singledetector \glsSNR of 10.
The corresponding error in for individual systems is consistent with that found in the available literature Damour et al. (2012); Lackey et al. (2012); Agathos et al. (2015). While these studies probe how accurately can be measured from \glsGW observations, they do not explore the effect of directly including inputs from nuclear theory, which is the point of this work.
To translate these correlated uncertainties in observables s and s (assuming effects of component spins to be small for ) to nuclear physics parameters, the \glsFIM generated from the waveforms described above is transformed to the space of nuclear parameters via the Jacobian such as the partial derivative . These are then combined with a \glsFIM from the base nuclear uncertainties, and information about neutron star masses and radii at levels expected of \glsNASA’s \glsNICER mission to obtain a final covariance matrix for the 18 parameters.
The Fisher method for estimating errors has limitations, one of the main being the need for a high \glsSNR. Bayesian methods are more reliable, but computationally much more expensive. For this latter reason use Fisher methods, whose computationally efficiency allows us to reduce source selection effects on the error estimates. We are able to quickly compute the \glsFIM for hundreds of binaries, characterizing the variance within the population models. In spite of the drawbacks, the Fisher errors quoted here make the case to invest in Bayesian methods.
Methodology
For a given population synthesis model, we simulate ten thousand \glsBNS systems and distribute them uniformly in comoving volume between a luminosity distance of and . The latter limit is not too far from the horizon distance () of the network of \glsaLIGO and Advanced Virgo detectors beyond which \glsBNS sources will produce signals with network \glsSNR of less than 8. Also, below we expect almost an order of magnitude fewer sources than those up to a distance of . This fact notwithstanding the measurement precision for a nearby source (\glsGW170817 was at a distance of ) can rival that of a population of more distant sources. This is why we also present results for a \glsGW170817like source at aLIGO design sensitivity.
Our main results are summarized in fig. 1, which shows how the constraints on the pressure of pure neutron matter improve as a function of additional \glsNICER or \glsLIGO observations. We start from the errors listed in table 1 which, for the purposes of this analysis, we interpret as uncorrelated normal errors for the parameters. In general, errors have been overestimated to ensure that our results are conservative. The resulting \glsFIM – a diagonal matrix of the inverse variances – provides our starting point. From this \glsFIM, we use forward error propagation to determine the error in pressure which we label “Nuclear”.
The largest uncertainty comes from the form of the \glsEoS in the core of the neutron star. Although a description in terms of homogeneous nuclear matter may persist to some depth, it is likely that there is some sort of phase transition to hyperonic or strange quark matter. The core \glsEoS is thus largely unknown. To assess the impact of large variations in the core \glsEoS, we compare the constraints obtained under a rather large variation of the core parameters, as well as in the presence of a strong firstorder phase transition (\glsCentral2c_trans). This comparison was summarized in fig. 5. Here we see rather large sensitivity to a smaller as expected: if this is small, the core transition occurs at low density, and not enough conventional nuclear matter exists to be sensitive to \glsGW observations. As long as the core transition is above or so, the constraints on are relatively insensitive to the form of the core \glsEoS unless there is a strong firstorder phase transition.
EoS
Vii Conclusion
The \glsGW170817 event demonstrated that useful constraints on the neutron star structure can be obtained from \glsplGW. In this article we have addressed how future observations can provide more detailed constraints on the properties of dense matter. By separating the neutron star into four distinct regions, and providing a unique nuclear physics based parameterization of the \glsEoS of the crust and outercore, we have analyzed how measurements of the tidal deformability can constrain nuclear properties of dense matter. Our parameterization, which uses the same underlying \glsEoS of neutron matter in both inner crust and outer core, allowed us to estimate for the first time constraints on the \glsEoS of pure neutron matter in the density interval where controlled calculations are becoming feasible. These constraints, as they become available, will provide valuable guidance for nuclear physics. In the inner core, where the \glsEoS is poorly constrained, the speed of sound is allowed to vary over a large range constrained only by causality and the requirement that the \glsEoS produce a 2 solar mass neutron star. We have taken first steps to study how the large uncertainties associated with the \glsEoS of the inner core limits our ability to constrain the \glsEoS of neutron matter in the outer core. The results we obtain suggest that, in the absence of strong firstorder transitions in the core, even a handful of detections can constrain the pressure of neutron matter in the density interval between and to better than 20%.
The principal component analysis presented in fig. 7 suggests that future \glsLIGO observations will provide strong constraints on the density dependence of the pure neutron matter \glsEoS in the outer core. In particular, we find that the exponent in the neutron matter \glsEoS defined in eq. 2 will be well constrained. As expected, the nuclear physics parameters are better constrained when the outer core makes the dominant contribution to the tidal deformability. This is the case when the neutron matter \glsEoS is stiff in the dense regions of the outer core and for lowmass neutron stars. If instead, the \glsEoS in the outer core is soft or if a strong firstorder phase transition were to occur at relatively lowdensity, constraints on the neutron matter \glsEoS are weaker. In these cases, the inner core has a larger impact on the tidal deformability and \glsGW detections will provide constraints for matter encountered in the inner core.
Although our focus here was to study the impact of the most common events that occur at large distances, we find that a single close by event similar to \glsGW170817 at at design sensitivity will provide valuable constraints. However, in the absence of such a nearby event, similar constraints may be realized by a dozen or so more distant events.
One limitation of our study is the simple parameterization of the \glsEoS of the inner core. While this is adequate as a first step, to constrain the \glsEoS of the inner core, a parameterization that allows for larger variability at high density will be needed. In addition, to gain more confidence in the constraints we have presented for neutron matter, it will be necessary to systematically marginalize over population models for neutron star masses and spins, and the uncertainty in the \glsEoS of the inner core. A Bayesian approach would be better suited for this purpose, and we are in the processes of developing computer programs needed for such a study.
Acknowledgements.
We thank K. G. Arun for helpful discussions and early collaboration on waveform models with tidal corrections. We also thank Philippe Landry for carefully reading the manuscript and making useful comments. \myscSb acknowledges partial support from the Navajbai Ratan Tata Trust. \myscSr acknowledges support from the \myscus Department of Energy Grant No. \myscdefg0200er41132 and from the National Science Foundation Grant No. \myscphy1430152 (\myscjina Center for the Evolution of the Elements). \myscAm acknowledges partial support from the \myscserb StartUp Research for Young Scientists Scheme project Grant No. \myscsb/ftp/ps067/2014, \myscdst, India.References
 Aasi et al. (2015) J. Aasi et al. (LIGO Scientific), Class. Quant. Grav. 32, 074001 (2015), arXiv:1411.4547 [grqc] .
 Acernese et al. (2015) F. Acernese et al. (VIRGO), Class. Quant. Grav. 32, 024001 (2015), arXiv:1408.3978 [grqc] .
 Abbott et al. (2017a) B. P. Abbott et al. (LIGO Scientific, Virgo), Phys. Rev. Lett. 119, 161101 (2017a), arXiv:1710.05832 [grqc] .
 Abbott et al. (2017b) B. P. Abbott et al. (LIGO Scientific, Virgo, Fermi GBM, INTEGRAL, IceCube, AstroSat Cadmium Zinc Telluride Imager Team, IPN, InsightHxmt, ANTARES, Swift, AGILE Team, 1M2H Team, Dark Energy Camera GWEM, DES, DLT40, GRAWITA, FermiLAT, ATCA, ASKAP, Las Cumbres Observatory Group, OzGrav, DWF (Deeper Wider Faster Program), AST3, CAASTRO, VINROUGE, MASTER, JGEM, GROWTH, JAGWAR, CaltechNRAO, TTUNRAO, NuSTAR, PanSTARRS, MAXI Team, TZAC Consortium, KU, Nordic Optical Telescope, ePESSTO, GROND, Texas Tech University, SALT Group, TOROS, BOOTES, MWA, CALET, IKIGW Followup, H.E.S.S., LOFAR, LWA, HAWC, Pierre Auger, ALMA, Euro VLBI Team, Pi of Sky, Chandra Team at McGill University, DFN, ATLAS Telescopes, High Time Resolution Universe Survey, RIMAS, RATIR, SKA South Africa/MeerKAT), Astrophys. J. 848, L12 (2017b), arXiv:1710.05833 [astroph.HE] .
 Abbott et al. (2017c) B. P. Abbott et al. (LIGO Scientific, Virgo, 1M2H, Dark Energy Camera GWE, DES, DLT40, Las Cumbres Observatory, VINROUGE, MASTER), Nature 551, 85 (2017c), arXiv:1710.05835 [astroph.CO] .
 Chen et al. (2018) H.Y. Chen, M. Fishbach, and D. E. Holz, Nature 562, 545 (2018), arXiv:1712.06531 [astroph.CO] .
 Nair et al. (2018) R. Nair, S. Bose, and T. D. Saini, Phys. Rev. D98, 023502 (2018), arXiv:1804.06085 [astroph.CO] .
 SoaresSantos et al. (2019) M. SoaresSantos et al. (DES, LIGO Scientific, Virgo), Submitted to: Astrophys. J. (2019), arXiv:1901.01540 [astroph.CO] .
 Abbott et al. (2018a) B. P. Abbott et al. (LIGO Scientific, Virgo), Phys. Rev. Lett. 121, 161101 (2018a), arXiv:1805.11581 [grqc] .
 Flanagan and Hinderer (2008) E. E. Flanagan and T. Hinderer, Phys. Rev. D77, 021502 (2008), arXiv:0709.1915 [astroph] .
 De et al. (2018) S. De, D. Finstad, J. M. Lattimer, D. A. Brown, E. Berger, and C. M. Biwer, (2018), arXiv:1804.08583 [astroph.HE] .
 Tews et al. (2018a) I. Tews, J. Margueron, and S. Reddy, (2018a), arXiv:1804.02783 [nuclth] .
 Abbott et al. (2019) B. P. Abbott et al. (LIGO Scientific, Virgo), Phys. Rev. X9, 011001 (2019), arXiv:1805.11579 [grqc] .
 Tews et al. (2018b) I. Tews, J. Carlson, S. Gandolfi, and S. Reddy, Astrophys. J. 860, 149 (2018b), arXiv:1801.01923 [nuclth] .
 Gandolfi et al. (2009) S. Gandolfi, A. Y. Illarionov, K. E. Schmidt, F. Pederiva, and S. Fantoni, Phys. Rev. C 79, 054005 (2009).
 Gandolfi et al. (2010) S. Gandolfi, A. Y. Illarionov, S. Fantoni, J. C. Miller, F. Pederiva, and K. E. Schmidt, MNRAS 404, L35 (2010).
 Gandolfi et al. (2012) S. Gandolfi, J. Carlson, and S. Reddy, Phys. Rev. C 85, 032801 (2012), arXiv:1101.1921 .
 Gandolfi et al. (2014a) S. Gandolfi, J. Carlson, S. Reddy, A. Steiner, and R. Wiringa, Eur. Phys. J. A 50, 10 (2014a), arXiv:1307.5815 [nuclth] .
 Miller and Lamb (2016) M. C. Miller and F. K. Lamb, Eur. J. Phys. A 52, 63 (2016).
 Dominik et al. (2012) M. Dominik, K. Belczynski, C. Fryer, D. Holz, E. Berti, T. Bulik, I. Mandel, and R. O’Shaughnessy, Astrophys. J. 759, 52 (2012), arXiv:1202.4901 [astroph.HE] .
 Miller (2016) M. C. Miller, Astrophys. J. 822, 27 (2016).
 Agathos et al. (2015) M. Agathos, J. Meidam, W. Del Pozzo, T. G. F. Li, M. Tompitak, J. Veitch, S. Vitale, and C. V. D. Broeck, Phys. Rev. D92, 023012 (2015), arXiv:1503.05405 [grqc] .
 (23) Supplementary Material.
 Horowitz et al. (2012) C. J. Horowitz, Z. Ahmed, C.M. Jen, A. Rakhman, P. A. Souder, M. M. Dalton, N. Liyanage, K. D. Paschke, K. Saenboonruang, R. Silwal, G. B. Franklin, M. Friend, B. Quinn, K. S. Kumar, D. McNulty, L. Mercado, S. Riordan, J. Wexler, R. W. Michaels, and G. M. Urciuoli, Phys. Rev. C 85, 032501 (2012).
 Horowitz et al. (2014a) C. Horowitz, K. Kumar, and R. Michaels, Eur. Phys. J. A 50, 48 (2014a), 10.1140/epja/i2014140483.
 Horowitz et al. (2014b) C. J. Horowitz, E. F. Brown, Y. Kim, W. G. Lynch, R. Michaels, A. Ono, J. Piekarewicz, M. B. Tsang, and H. H. Wolter, J. Phys. G 41, 093001 (2014b).
 Hebeler and Schwenk (2010) K. Hebeler and A. Schwenk, Phys. Rev. C 82, 014314 (2010).
 Hebeler et al. (2013) K. Hebeler, J. M. Lattimer, C. J. Pethick, and A. Schwenk, Astrophys. J. 773, 11 (2013).
 Wlazłowski et al. (2014) G. Wlazłowski, J. W. Holt, S. Moroz, A. Bulgac, and K. J. Roche, Phys. Rev. Lett. 113, 182503 (2014), arXiv:1403.3753 .
 Gandolfi et al. (2014b) S. Gandolfi, A. Lovato, J. Carlson, and K. E. Schmidt, Phys. Rev. C 90, 061306 (2014b), arXiv:1406.3388 [nuclth] .
 Lynn et al. (2015) J. E. Lynn, I. Tews, J. Carlson, S. Gandolfi, A. Gezerlis, K. E. Schmidt, and A. Schwenk, “Chiral threenucleon interactions in light nuclei, neutron scattering, and neutron matter,” (2015), arXiv:1509.03470 .
 Page and Reddy (2006) D. Page and S. Reddy, Annu. Rev. Nucl. Part. Sci. 56, 327 (2006), arXiv:astroph/0608360v1 .
 Steiner and Gandolfi (2012) A. W. Steiner and S. Gandolfi, Phys. Rev. Lett. 108, 081102 (2012).
 Steiner et al. (2013) A. W. Steiner, J. M. Lattimer, and E. F. Brown, Astrophys. J. Lett. 765, L5 (2013).
 Lattimer and Steiner (2014) J. M. Lattimer and A. W. Steiner, Eur. Phys. J. A 50, 1 (2014).
 Baym et al. (1971) G. Baym, C. Pethick, and P. Sutherland, Astrophys. J. 170, 299 (1971).
 Negele and Vautherin (1973) J. W. Negele and D. Vautherin, Nucl. Phys. A 207, 298 (1973).
 Sharma et al. (2015) B. K. Sharma, M. Centelles, X. Vinas, M. Baldo, and G. F. Burgio, Astron. & Astrophys. 584, A103 (2015), arXiv:1506.00375 .
 Haensel et al. (2007) P. Haensel, A. Y. Potekhin, and D. G. Yakovlev, Neutron Stars 1, 1st ed., Astrophysics and Space Science Library, Vol. 326 (SpringerVerlag, New York, 2007).
 Chamel and Haensel (2008) N. Chamel and P. Haensel, Living Rev. Relativity 11 (2008), 10.12942/lrr200810, arXiv:0812.3955 .
 Nelson et al. (2018) A. Nelson, S. Reddy, and D. Zhou, (2018), arXiv:1803.03266 [hepph] .
 Read et al. (2009) J. S. Read, B. D. Lackey, B. J. Owen, and J. L. Friedman, Phys. Rev. D 79, 124032 (2009).
 Fortin et al. (2016) M. Fortin, C. Providência, A. R. Raduta, F. Gulminelli, J. L. Zdunik, P. Haensel, and M. Bejger, Phys. Rev. C 94, 035804 (2016), arXiv:1604.01944 .
 Zdunik et al. (2016) J. L. Zdunik, M. Fortin, and P. Haensel, “Neutron star properties and the equation of state for its core,” (2016), arXiv:1611.01357 .
 Ravenhall et al. (1983) D. G. Ravenhall, C. J. Pethick, and J. R. Wilson, Phys. Rev. Lett. 50, 2066 (1983).
 Chamel et al. (2007) N. Chamel, S. Naimi, E. Khan, and J. Margueron, Phys. Rev. C 75, 055806 (2007).
 Lattimer et al. (1985) J. Lattimer, C. Pethick, D. Ravenhall, and D. Lamb, Nucl. Phys. A 432, 646 (1985).
 Steiner (2012) A. W. Steiner, Phys. Rev. C 85, 055804 (2012).
 Margueron et al. (2018) J. Margueron, R. Hoffmann Casali, and F. Gulminelli, Phys. Rev. C 97, 025805 (2018).
 Gandolfi et al. (2015) S. Gandolfi, A. Gezerlis, and J. Carlson, Annu. Rev. Nucl. Part. Sci. 65, 303 (2015), arXiv:1501.0567 .
 Roggero et al. (2014) A. Roggero, A. Mukherjee, and F. Pederiva, Phys. Rev. Lett. 112, 221103 (2014).
 Rrapaj et al. (2016) E. Rrapaj, A. Roggero, and J. W. Holt, prc 93, 065801 (2016).
 Lee et al. (1998) C.H. Lee, T. T. S. Kuo, G. Q. Li, and G. E. Brown, Phys. Rev. C 57, 3488 (1998).
 GonzalezBoquera et al. (2017) C. GonzalezBoquera, M. Centelles, X. Viñas, and A. Rios, Phys. Rev. C 96, 065806 (2017).
 Li et al. (2008) B.A. Li, L.W. Chen, and C. M. Ko, Phys. Rep. 464, 113 (2008).
 Bulgac et al. (2018) A. Bulgac, M. M. Forbes, S. Jin, R. N. Perez, and N. Schunck, Phys. Rev. C 97, 044313 (2018), arXiv:1708.08771 [nuclth] .
 (57) M. Alford, Raised in discussions at the INT162b program.
 Vines et al. (2011) J. Vines, E. E. Flanagan, and T. Hinderer, Phys. Rev. D83, 084051 (2011), arXiv:1101.1673 [grqc] .
 Del Pozzo et al. (2013) W. Del Pozzo, T. G. F. Li, M. Agathos, C. Van Den Broeck, and S. Vitale, Physical Review Letters 111, 071101 (2013), arXiv:1307.8338 [grqc] .
 Bose et al. (2018) S. Bose, K. Chakravarti, L. Rezzolla, B. S. Sathyaprakash, and K. Takami, Phys. Rev. Lett. 120, 031102 (2018), arXiv:1705.10850 [grqc] .
 Abbott et al. (2018b) B. P. Abbott et al. (LIGO Scientific, Virgo), (2018b), arXiv:1811.12907 [astroph.HE] .
 Stovall et al. (2018) K. Stovall et al., Astrophys. J. 854, L22 (2018), arXiv:1802.01707 [astroph.HE] .
 Buonanno et al. (2009) A. Buonanno, B. Iyer, E. Ochsner, Y. Pan, and B. S. Sathyaprakash, Phys. Rev. D80, 084043 (2009), arXiv:0907.0700 [grqc] .
 Dietrich et al. (2017) T. Dietrich, S. Bernuzzi, and W. Tichy, Phys. Rev. D96, 121501 (2017), arXiv:1706.02969 [grqc] .
 Hinderer (2008) T. Hinderer, Astrophys. J. 677, 1216 (2008), arXiv:0711.2420 [astroph] .
 Postnikov et al. (2010) S. Postnikov, M. Prakash, and J. M. Lattimer, Phys. Rev. D 82, 024016 (2010), arXiv:1004.5098 .
 Ajith and Bose (2009) P. Ajith and S. Bose, Phys. Rev. D79, 084032 (2009), arXiv:0901.4936 [grqc] .
 Helstrom (1995) C. W. Helstrom, Elements of signal detection and estimation (PrenticeHall, Inc., Upper Saddle River, NJ, USA, 1995).
 aLI (2010) Advanced LIGO anticipated sensitivity curves, Tech. Rep. LIGOT0900288v3 (LIGO Scientific Collaboration, https://dcc.ligo.org/LIGOT0900288/public, 2010).
 Vallisneri (2008) M. Vallisneri, Phys. Rev. D 77, 042001 (2008), arXiv:grqc/0703086 .
 Rodriguez et al. (2013) C. L. Rodriguez, B. Farr, W. M. Farr, and I. Mandel, Phys. Rev. D88, 084013 (2013), arXiv:1308.1397 [astroph.IM] .
 Rodriguez et al. (2014) C. L. Rodriguez, B. Farr, V. Raymond, W. M. Farr, T. B. Littenberg, D. Fazi, and V. Kalogera, Astrophys. J. 784, 119 (2014), arXiv:1309.3273 [astroph.HE] .
 Damour et al. (2012) T. Damour, A. Nagar, and L. Villain, Phys. Rev. D85, 123007 (2012), arXiv:1203.4352 [grqc] .
 Lackey et al. (2012) B. D. Lackey, K. Kyutoku, M. Shibata, P. R. Brady, and J. L. Friedman, Phys. Rev. D85, 044061 (2012), arXiv:1109.3402 [astroph.HE] .
Appendix A Supplementary Material
a.1 Surface Term in the CLDM
In our implementation of the \glsCLDM, we use the following surface term with an effective surface tension following Lattimer et al. (1985) (see also Steiner (2012)): \cref@addtoresetequationparentequation
(9a)  
(9b)  
where is the proton fraction in the nucleus, and . We parameterize this as where is held fixed as a parameter of the theory, and is varied to smoothly match the tabulated outercrust data. 
a.2 Polytropes
In fig. 10 we compare the constraints obtained on the total pressure of nuclear matter in equilibrium using our \glsCentral2c unified parameterization with those obtained using the piecewise polytropic \glsEoS in Read et al. (2009). To better compare these, we do the following:

Fit the parameters of the polytrope to best match our \glsCentral2c \glsEoS: , , , and .

We use the same speedofsound core parameterization with , , as our \glsCentral2c \glsEoS.

We start with a bare “Nuclear” constraint by computing the covariance matrix of the parameters from Table III of Read et al. (2009) over the following \glsEoS models that have a small pressure at saturation density: \myscpal6, \myscsly, \myscapr1, \myscapr2, \myscapr3, \myscapr4, \myscfps, \myscwff1, \myscwff2, \myscwff3, \myscbbb2, \myscbpal12, \mysceng, \myscmpa1, \myscbgn1h1, \myscpcl2, \myscalf1, \myscalf2, \myscalf3, and \myscalf4. (This excludes some models with hyperon (\myscgnh3, \mysch17), pion (\myscps), and kaon (\myscgs12) condensates, as well as the strangequark matter models \myscms12, which all have significantly higher saturation pressures ). This gives similar bare “Nuclear” errors as our \glsCentral2c model at and above saturation density.
We note that the constraints on are very similar to those from our “Nuclear” parameter set. To obtain this, however, it was critical to use correlated errors in the polytrope parameters. To this end, taking a polytropic \glsEoS with uncorrelated priors is inadvisable. Only once correlated priors are used does the polytropic equation of state provide constraints comparable to those that can be obtained from our nuclear parameterization.
a.3 Tabulated EoS Data
Here we present somewhat tighter constraints on tabulated \glsEoS data, required to ensure convexity, than we have seen presented in the literature. Suppose we have an interval with tabulated density, pressure, and energy , , and . If these data come from an equation of state that satisfies thermodynamic convexity and causality , then each interval must satisfy the following conditions:
(10) 
The tabulated date in Sharma et al. (2015) used for the outer crust required some minor corrections to ensure these constraints are met.
a.4 Thermodynamic Relationships
Here we briefly review some thermodynamic relationships for an \glsEoS with a single conserved component with density and chemical potential , energy density , energy per particle , and pressure . These are used at various places throughout the text, such as relating the slope of the symmetry energy to the pressure of neutron matter in eq. 1b \cref@addtoresetequationparentequation
(11a)  
(11b)  
(11c) 
a.5 Comparison Plots
On the following pages, we provide comparison plots for all of the \glsEoS models discussed in the text.