Constraining the neutron-matter equation of state with gravitational waves

Constraining the neutron-matter equation of state with gravitational waves

Michael McNeil Forbes Department of Physics & Astronomy, Washington State University, Pullman, Washington 99164–2814, \myscusa Department of Physics, University of Washington, Seattle, Washington 98195–1560, \myscusa    Sukanta Bose Department of Physics & Astronomy, Washington State University, Pullman, Washington 99164–2814, \myscusa Inter-University Centre for Astronomy and Astrophysics, Post Bag 4, Ganeshkhind, Pune 411 007, India    Sanjay Reddy Institute for Nuclear Theory, University of Washington, Seattle, Washington 98195–1560, \myscusa    Dake Zhou Department of Physics, University of Washington, Seattle, Washington 98195–1560, \myscusa    Arunava Mukherjee Max-Planck-Institut für Gravitationsphysik (Albert Einstein Institute), D-30167 Hannover, Germany Leibniz Universität Hannover, D-30167 Hannover, Germany    Soumi De Department of Physics, Syracuse University, Syracuse, New York 13244, \myscusa

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 neutron-matter \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 x-ray 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 so-called parabolic approximation to extrapolate from symmetric nuclear matter, allowing us to directly constrain the neutron-matter \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 population-model sensitivity, phase transitions in the core \glsEoS, and large deviations from the central parameter values.

preprint: \myscligo-p1900097preprint: \myscint-pub-19-009

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 multi-messenger 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); Soares-Santos 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 many-body 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 low-density and symmetric properties of nuclear matter to significantly improve constraints on neutron-rich matter at higher densities.

In this article we discuss how we can extract more detailed information about the properties of dense neutron-rich matter and neutron stars during the next few years with more \glsGW detections and measurements of neutron star radii expected from x-ray 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 many-body nuclear forces. This, combined with advanced computational methods to solve the non-relativistic quantum many-body 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 neutron-rich 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.

Figure 1: \glsresetEoS(color online) Relative constraints on the pressure of neutron matter from simulated merger events, and expected constraints from \glsNICER (J0437) Miller and Lamb (2016) (, ). From top: constraints from nuclear theory augmented by \glsNICER, from a single merger event at with \glsaLIGO sensitivity, then various \glsLIGO events drawn from \glsCompleteStdAsubsolarNSNS that have . The shading shows the range of sampling errors ( or 68th percentile) demonstrating variation within the \glsCompleteStdAsubsolarNSNS population model Dominik et al. (2012). Beyond the vertical yellow line, we use the core \glsEoS. Inset: with error bands corresponding to each of the constraints.

Ii Results

Our main result is that even a handful of \glsGW observations of \glsBNS mergers will provide the most stringent constraints on the low-temperature 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.

  • An uncorrelated mass and radius measurement of \myscj0437 projected to be measured at a 5% level from \glsNASA’s \glsNICER mission – i.e. with a 10% measurement of  Miller (2016); Miller and Lamb (2016).

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 non-linear effects, we provide similar plots for comparison in EPA () for the different central values listed in table 2.

Figure 2: (color online) Improvement in relative constraint on the pressure of neutron matter at and (related to the slope of the symmetry energy) to an increasing number of simulated merger events applied to the initial nuclear constraints denoted with a plus at . The shading shows the range of sampling errors ( or 68th percentile) demonstrating variation within the \glsCompleteStdAsubsolarNSNS population model. The lower dotted curve shows the level of the most tightly constrained principal component (1st \myscpc).

To put these results in perspective, consider the nuclear symmetry energy and the slope of its density dependence , \cref@addtoresetequationparentequation


where is the energy-per-particle of uniform nuclear matter. If the so-called 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 .



Figure 3: (color online) Regions of the neutron star. The upper three wedges represent a cross-section of , , and neutron stars respectively. As discussed in the text, the unified \acrshortEoS smoothly connects four distinct regions from low density on the left to high density on the right. The radius of these transitions for the \glsCentral2c parameter values is shown in the top plot. These are connected to the equation of state expressed in terms of the pressure (solid (black) line on left axis) as a function of the total baryon density in units of the saturation density . From low to high density, the regions of the \acrshortEoS are: a) the outer crust (very low density which is too small to see on the lower plot) that interpolates the data of Baym et al. (1971) and Negele and Vautherin (1973) as tabulated in Sharma et al. (2015) (blue) with minor corrections to ensure convexity as discussed in EPA (); b) the inner crust modeled by the \acrshortCLDM Haensel et al. (2007); *Chamel:2008 (orange); c) the outer core of homogeneous nuclear matter in beta-equilibrium (green); d) the inner core equation of state parameterized by a quadratic speed of sound (red). At the right, the various (red) dashed lines correspond to the core density of the respective stars. At the bottom are corresponding dashed curves (purple) proportional to (normalized to the maximum value on the right axis) for the two lower-mass stars. This roughly correlates with the local contribution to the dimensionless tidal deformability Nelson et al. (2018).

Iii Parameterization of the Nuclear Equation of State


EoS 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 neutron-rich 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 neutron-rich 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 speed-of-sound parameterization discussed in section III.3.

Figure 4: (color online) Mass-radius curves for the \glsplEoS considered in table 2. The thick solid curve is our \glsCentral2c \glsEoS. Dashed curves correspond to different core parameterizations. Thin curves correspond to \glsplEoS for which astrophysical observations would provide poor constraints for nuclear physics. These include a sharp first-order transition in the core (\glsCentral2c_trans), and soft \glsplEoS (\glsSoft2c and \glsStiff2c_Soft) which form very compact objects with low deformability. The \glsCentral2c_low_Ec \glsEoS also poorly constrains nuclear physics since the core appears close to the saturation density. As shown later in fig. 9, for these types of \glsEoS, observations constrain the core parameters rather than the properties of neutron matter.

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 low-mass 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


CLDM 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 Wigner-Seitz 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 Wigner-Seitz cell precludes the possibility of pasta phases Ravenhall et al. (1983) the errors incurred by the Wigner-Seitz 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 first-order 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:
Symmetric nuclear matter and symmetry parameters:
Neutron matter parameters:
Proton polaron parameters:
Inner-core parameters:
Table 1: Parameters defining the \glsCentral2c \glsEoS along with their uncorrelated covariance (expressed using the \glsSI convention ) used to defined the “Nuclear” error estimates prior to including information from astrophysical observations. We take the values of the \glsCLDM parameters from the fits to the \myscapr \glsEoS tabulated in Steiner (2012) but assign large errors to encompass missing physics such as the possibility of pasta phases. Symmetric nuclear matter and symmetry parameters have errors taken from the extensive analysis Margueron et al. (2018). Neutron matter parameters have errors estimated from \glsQMC calculations with various three-body interactions Gandolfi et al. (2009); *Gandolfi:2010b; *Gandolfi:2012; *Gandolfi:2014a, and are consistent with recent \glsQMC results based on chiral \glsEFT interactions Wlazłowski et al. (2014); Gandolfi et al. (2014b, 2015); Lynn et al. (2015). Proton polaron parameters have errors estimated from the \glsQMC calculations Roggero et al. (2014) and are consistent with estimates from chiral interactions Rrapaj et al. (2016). The core parameters are chosen to allow for a star at the extremes of all of our models except for the \glsSoft2c \glsEoS which requires a lower core transition and are given large errors to be conservative with the exception of the parameter . This is given a small error for the purposes of our statistical analysis as the dependence is highly non-linear. Variations of this parameter are considered specially in fig. 5.
Neutron Matter Inner Core
\glsEoS [] [] []
Table 2: List of changed \glsEoS parameters compared in this work. All other parameters share the same values as the \glsCentral2c \glsEoS in the top row, which takes the central values listed in table 1. The first four variations – \glsSoft2c, \glsStiff2c, \glsSoft2c_Stiff, and \glsStiff2c_Soft – refer to the properties of the neutron-matter equation of state and whether the \glsEoS of the outer core is softer or stiffer than \glsCentral2c at low/high density. The next three variations – \glsCentral2c_low_Ec, \glsCentral2c_high_Ec, and \glsCentral2c_low_C_max – explore variations of the core \glsEoS. 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 first-order phase transition with discontinuity . (See fig. 5.)
Figure 5: (color online) Sensitivity of the constraint on the pressure of neutron matter from simulated merger events drawn from the \glsCompleteStdAsubsolarNSNS population model to variations of the core equation of state. The vertical yellow lines denote the density at which the \glsEoS reverts to the core form. Inset: form of the various core speed-of-sound functions .

iii.2 Homogeneous Nuclear Matter

One of the main new features of our analysis is to parameterize the nuclear-matter \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 so-called parabolic approximation, which is valid only if quadratic terms dominate over quartic and higher-order terms. While there is some support for this below saturation density from relativistic \glsDBHF calculations Lee et al. (1998), Gogny forces Gonzalez-Boquera 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 neutron-matter \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:


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 three-body 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:


where is the proton effective mass, and describes the self-energy of the proton polaron. This function is presently poorly constrained by \glsQMC and experimental data and all known results are consistent with a simple two-parameter quadratic expansion:


where and where returns to zero. (We expect to curve up for higher densities due to the repulsive nature of nuclear three-body 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 :


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 high-density \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 first-order 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 point-mass 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, tens-to-hundreds 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 non-spinning 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


where and denote its amplitude and phase in the time domain. For \glsFIM-based parameter estimation, we work with the Fourier transform of the strain above. This is constructed by adding to the point-particle 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 cross-section 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 sub-solar metallicity. Component masses are narrowly peaked for solar metallicity systems while subsolar metalicity system produce a wider mass distribution.

Figure 6: \glsresetEoS(color online) Population model sensitivity of the constraint on the pressure of neutron matter from simulated merger events drawn from various different population models. The weaker constraints from the \glsUniform model result from distributing the events over larger mass objects. As shown in fig. 8, this provides more information about the properties of the core at the expense of information about the lower-density regions that constrain neutron matter.

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


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


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


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


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


Uniform sampling of neutron stars with masses between and .

Figure 7: (color online) Principal component analysis of the simulated observational data in terms of the \glsEoS parameters. Each column is a plot of the components of most significantly constrained eigenvector for the particular combination of observations listed at the bottom. These should be interpreted as follows: A linear combination of the log of the corresponding parameters is constrained to the tolerance shown at the top. The rightmost column shows the principal component analysis for simulated merger events drawn from the \GlsCompleteStdAsubsolarNSNS population model, and is the same as the leftmost column of fig. 9. The errors in the tolerances, shown as small black strips in middle of the component bars, are obtained by performing 200 independent samples and demonstrate variation within the population model. (These errors are small here, but quite visible in the second principal components of fig. 8.)
Figure 8: (color online) First two principal components for \glsGW observations drawn from the \GlsCompleteStdAsubsolarNSNS (left) and \GlsUniform (right) population models. This demonstrates the wider distribution of masses in the \glsUniform model as compared to \GlsCompleteStdAsubsolarNSNS. The narrow distribution in \GlsCompleteStdAsubsolarNSNS leads to tighter statistical constraints on the 1st principal component, but leaves other directions in parameter space poorly explored. In contrast, the \glsUniform model distributes the 15 events over a larger range of masses, reducing the constraints on the 1st principle component, but providing more information about a second direction. Even for observations, the next principal is poorly constrained at a level worse than 200%: more observations would be required to constrain this component at a useful level. Thus, neutron-star observables seem to provide tight constraints in a single direction of parameter space.
Figure 9: (color online) Principal component analysis of simulated merger events drawn from the \glsCompleteStdAsubsolarNSNS population model for each of the \glsEoS parameters listed in table 2. The leftmost column thus corresponds to the rightmost column of fig. 7. This analysis makes clear the non-linear dependence of the problem on the \glsEoS parameters: neutron star observations constrain either properties of the core for parameter values such as \glsSoft2c, \glsStiff2c_Soft, \glsCentral2c_low_Ec, or \glsCentral2c_trans, or the neutron matter \glsEoS for more central values.

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 frequency-domain 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 matched-filter 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 noise-limited 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 variance-covariance matrix,  Helstrom (1995).

The matrix is bounded by the signal via the Cramer-Rao inequality, which states that


where is the \glsFIM:


Above, is the partial derivative with respect to the parameter and is the one-sided noise \glsPSD Helstrom (1995). We take the latter to be the zero-detuned high-power \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 error-estimates for total-mass and mass-ratio (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 single-detector \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.


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 \glsGW170817-like 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 over-estimated 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 first-order 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 first-order phase transition.



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 outer-core, 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 first-order 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 low-mass neutron stars. If instead, the \glsEoS in the outer core is soft or if a strong first-order phase transition were to occur at relatively low-density, 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.

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. \myscde-fg02-00er41132 and from the National Science Foundation Grant No. \myscphy-1430152 (\myscjina Center for the Evolution of the Elements). \myscAm acknowledges partial support from the \myscserb Start-Up Research for Young Scientists Scheme project Grant No. \myscsb/ftp/ps-067/2014, \myscdst, India.


Figure 10: (color online) Relative constraints on the pressure of nuclear matter in -equilibrium for the \glsCentral2c \glsEoS (dark curves) and a polytropic \glsEoS (light curves) with the same form as Read et al. (2009) with parameters fit to give a similar mass-radius curve, but using the same core \glsEoS as ours. (Note: Kinks in these curves occur when the form of the \glsEoS changes - for example, just below , the from of the polytrope Read et al. (2009) changes.)

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

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 outer-crust 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:

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

  2. We use the same speed-of-sound core parameterization with , , as our \glsCentral2c \glsEoS.

  3. 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, \mysch1-7), pion (\myscps), and kaon (\myscgs1-2) condensates, as well as the strange-quark matter models \myscms1-2, 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:


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


a.5 Comparison Plots

On the following pages, we provide comparison plots for all of the \glsEoS models discussed in the text.

Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description