Reconstructing the neutronstar equation of state with gravitationalwave detectors from a realistic population of inspiralling binary neutron stars
Abstract
Gravitationalwave observations of inspiralling binary neutron star systems can be used to measure the neutronstar equation of state (EOS) through the tidally induced shift in the waveform phase that depends on the tidal deformability parameter . Previous work has shown that , a function of the neutronstar EOS and mass, is measurable by Advanced LIGO for a single event when including tidal information up to the merger frequency. In this work, we describe a method for stacking measurements of from multiple inspiral events to measure the EOS. We use Markov chain Monte Carlo simulations to estimate the parameters of a 4parameter piecewise polytrope EOS that matches theoretical EOS models to a few percent. We find that, for “realistic” event rates ( binary neutron star inspiral events per year with signaltonoise ratio in a single Advanced LIGO detector), combining a year of gravitationalwave data from a threedetector network with the constraints from causality and recent high mass neutronstar measurements, the EOS above nuclear density can be measured to better than a factor of two in pressure in most cases. We also find that in the mass range –, the neutronstar radius can be measured to better than km and the tidal deformability can be measured to better than g cm s (10%–50% depending on the EOS and mass). The overwhelming majority of this information comes from the loudest events. Current uncertainties in the postNewtonian waveform model, however, lead to systematic errors in the EOS measurement that are as large as the statistical errors, and more accurate waveform models are needed to minimize this error.
pacs:
97.60.Jd 26.60.Kp, 04.30.Tv,I Introduction
Observations of neutron stars (NSs), consisting of nuclear matter in the ground state with densities up to several times nuclear saturation density ( g/cm), in principle provide an ideal way to measure the nuclear equation of state (EOS) that describes the pressure as a function of density . At these densities, the most rigorous EOS constraints from NSs have come from measurements of high mass pulsars in compact binary systems, and it is now clear that the EOS must allow for NS masses Demorest et al. (2010); Antoniadis et al. (2013). More precise information can be obtained if the NS mass and radius are simultaneously measured, and Lindblom demonstrated explicitly that there is a onetoone map between the relations and Lindblom (1992).
However, various calculations of the mass and radius of NSs from available observations have produced only marginally consistent results. For example, Özel et al. used observations of thermonuclear bursts from xray binaries to measure the NS mass and radius, and found, for the three systems considered, 95% confidence intervals that were all in the range 9–12 km at (from Fig. 1 of Ref. Özel et al. (2010) before stacking observations). Steiner et al. also used mass and radius measurements from three xray bursters as well as from three quiescent lowmass xray binaries and found the NS radius to be 10.7–12.5 km at (from the range of 95% confidence intervals in Tables 7 and 8 of Ref Steiner et al. (2010)). Finally, Guillot et al. measured the mass and radius for five quiescent lowmass xray binaries, and found that the NS radius was 7.6–10.4km (90% confidence) assuming the radius was an approximately constant function of mass Guillot et al. (2013).
Gravitationalwave (GW) observations of NSs in inspiralling compact binaries, on the other hand, are sensitive only to the density profile of NSs, and are therefore free of the modeldependent uncertainties in the emission and absorption mechanisms that plague electromagnetic observations. In the next few years, observing runs will begin for second generation detectors, including the two Advanced LIGO (aLIGO) detectors Harry and LIGO Scientific Collaboration (2010) and Advanced Virgo (aVirgo) F. Acernese, et al. (Virgo Collaboration) (2009), and design sensitivity will likely be reached by the end of the decade LIGO Scientific Collaboration et al. (2013). In addition, KAGRA Somiya (2012) (formerly LCGT) and possibly LIGOIndia Iyer et al. (2011) will come online a few years later.
The EOS information provided by these detectors comes mainly from tidal interactions during the inspiral of binary neutron star (BNS) and black holeneutron star (BHNS) systems that induce quadrupolar deformations in the NSs Flanagan and Hinderer (2008). In the quasistationary approximation, the quadrupole moment of one star depends on the tidal field from the monopole of the other star through the relation . Here, is the EOS dependent tidal deformability and is related to the NS’s dimensionless Love number , first calculated in Ref. Hinderer (2008), and radius through the relation , where is the gravitational constant.
Because the tidal effect is a strong function of the mass ratio, and observed BNS events are likely to be significantly closer, this tidal interaction is likely to be seen for only BNS systems Pannarale et al. (2011). Additional information may also come from modes in the postmerger remnant of BNS systems Stergioulas et al. (2011); Clark et al. (2014) as well as the damping of quasinormal modes for low mass, high spin BHNS systems Lackey et al. (2014).
The measurability of tidal parameters in BNS systems with aLIGO for GW frequencies below 450Hz (prior to the last GW cycles before merger where higher order corrections to the tidal effect become important) was first examined for polytropic EOSs in Ref. Flanagan and Hinderer (2008) and for theoretical hadronic and quark matter EOSs in Ref. Hinderer et al. (2010). These studies found that for a single aLIGO detector, tidal interactions were only observable during this early inspiral stage for stiff EOSs, NS masses below 1.4 , and for rare, nearby sources with signaltonoise ratios (SNR) . Damour et al., however, found that when also including EOS information from the last GW cycles prior to contact, tidal parameters are in fact observable even with more common SNRs of Damour et al. (2012).
These studies relied on the Fisher matrix approximation for parameter estimation which assumes the distribution of parameters follows a multivariate Gaussian. Recent works using a fully Bayesian analysis for the threedetector aLIGOaVirgo network have confirmed that the tidal signal is indeed measurable when including the entire inspiral up to merger Del Pozzo et al. (2013); Wade et al. (2014). In particular Del Pozzo et al. used a method to stack tens of observations of a large number of BNS inspiral events to measure Del Pozzo et al. (2013). By parametrizing with a linear fit, they found that could be measured to at , but the mass dependence of could not be found. This was true even though they used an unrealistically large NS mass range of – for their simulated population.
In this paper we will demonstrate the advantages of parameterizing the EOS instead of . The main benefit comes from including prior knowledge of the EOS that would be more difficult to incorporate into the fit. As will be discussed in Section III, the requirements that the EOS is monotonic and causal provide the simple constraint that the slope of the EOS is in the range , where is the energy density and is the speed of light. Additionally, the constraints from the observed NSs can be included by simply rejecting any EOS parameters that lead to a maximum NS mass below . We will find that incorporating this information allows us to make significantly stronger statements about the EOS, radius, and tidal deformability, and this is true for a narrower, more realistic range of NS masses than simulated by Del Pozzo et al. Del Pozzo et al. (2013).
The BNS inspiral signal is almost exactly characterized by the mass, spin, and tidal deformability of each NS, and is therefore free from the intrinsic variability and uncertainties that bias electromagnetic measurements of the mass and radius. Currently unknown terms in the postNewtonian (PN) waveform model, however, will lead to systematic errors in recovering the tidal deformabilities. Using the Fisher matrix approximation Favata (2014); Yagi and Yunes (2014) as well as a Markov chain Monte Carlo (MCMC) Bayesian analysis Wade et al. (2014), it was found that the systematic error in from waveform uncertainties can be as large as the statistical error for a single BNS inspiral observation. When stacking observations the problem becomes worse, as the statistical errors decrease with more observations but the systematic errors do not, and we will discuss how systematic errors impact the recovery of the EOS.
We organize the paper as follows. We describe the BNS waveform model in Section II and the EOS parameterization in Section III. Our Bayesian method for estimating EOS parameters is derived in Section IV. We then show our results for a range of observation scenarios in Section V and the impact from waveform uncertainties in Section VI. Finally, we summarize our results in Section VII.
Ii The BNS inspiral signal
The strain observed in a detector from a GW is related to the two polarizations of the GW, and , by the detector’s antenna beam pattern response, and , through the relation
(1) 
and depend on the sky position (given by the right ascension and declination ) and GW polarization angle of the source.
For a BNS system in a quasicircular inspiral with an inclination angle , comoving (transverse) distance , and component masses and , the two polarizations of the GW are
(2)  
(3) 
where we have restricted the amplitude to the leading order and only keep the leading harmonic. Here, is the total mass, is the symmetric mass ratio, is the orbital phase, and is the standard PN order parameter defined by
The time evolution of and are evaluated from the PN expressions for the energy and luminosity via the energy balance requirement :
(4)  
(5) 
The energy and luminosity have the schematic form
(6)  
(7) 
where we write the tidal terms as functions of the dimensionless quantities .
The pointparticle PN terms for nonspinning systems depend only on and . The energy was recently calculated to 4PN order Bini and Damour (2013); Foffa and Sturani (2014), and the luminosity is known to 3.5PN order Blanchet (2014). For consistency with previous works Wade et al. (2014); Del Pozzo et al. (2013); Favata (2014); Yagi and Yunes (2014), we keep the pointparticle PN corrections to 3.5PN order.
The leading tidal terms begin at the same order as 5PN pointparticle terms, and although the tidal potential is known to the nexttonexttoleadingorder term Bini et al. (2012), we only include the leading and nexttoleading order terms for consistency with previous works Wade et al. (2014); Del Pozzo et al. (2013); Favata (2014); Yagi and Yunes (2014). These terms are Vines et al. (2011)
(8)  
(9) 
where .
As described in detail in Ref. Wade et al. (2014), there are several methods, first cataloged in Ref. Damour et al. (2001), to determine the phase evolution of the binary from the above expressions for its energy and luminosity, and we will use three of them. In the TaylorT1 method, one numerically integrates Eqs. (4) and (5) to find and up to the time and phase constants and . The TaylorT4 method is a slight variation where one reexpands the ratio in Eq. (5) then truncates the series at the highest known PN order (here 3.5PN) before numerically integrating. Finally in the TaylorF2 method, one approximates the Fourier transform of the waveform
(10) 
with the stationary phase approximation. The result, to leading order in the amplitude, is
(11)  
(12) 
Here, is the chirp mass, and, in the TaylorF2 approach, the phase has the analytic form
(13) 
where the pointparticle terms are provided in Ref. Buonanno et al. (2009). As was demonstrated in Refs. Flanagan and Hinderer (2008); Favata (2014); Wade et al. (2014), the individual tidal parameters and are highly correlated, so it is instead easier to reparameterize the tidal contribution to the phase in terms of the linear combinations
(14)  
(15) 
The tidal contribution then takes the simple form
(16) 
The two terms containing are significantly larger than the term containing , and, as was previously found, is not measurable with aLIGO Wade et al. (2014). This waveform model can therefore be expressed in terms of the 11 parameters .
For all versions of the PN waveform, we cut off the inspiral at the GW frequency corresponding to the Schwarzschild innermost stable circular orbit (ISCO) . However, for large NS radii corresponding to a stiff EOS, NSs can merge before reaching . For simplicity, we also choose to examine nonspinning BNS systems, although it has been shown that not including spin parameters can noticeably bias parameter estimation even for systems with relatively small spin magnitudes Favata (2014). A parallel investigation which studies how the choice of highfrequency cutoffs and spin can affect the estimation of tidal parameters is nearing completion Agathos and et al. (2014).
Iii The EOS
Because the EOS can be calculated from the Lindblom (1992) or relations Lindblom and Indik (2012, 2014), we can choose to parameterize either , or if we want to reconstruct the EOS. As Del Pozzo et al. found by parameterizing with a linear fit, it is difficult to accurately measure over a wide range of masses with BNS inspiral observations Del Pozzo et al. (2013). Instead, they found that could only be accurately measured for a specific fiducial mass which they chose to be . This is because, when parameterizing the function , there is little a priori information about the allowed functional form of . Lattimer and Steiner found similar results for mass and radius observations (Section 4.1 of Ref. Lattimer and Steiner (2014)) when they parameterized instead of the EOS. In contrast, it is much simpler to place useful constraints on the functional form of the EOS fit, and this allows one to make significantly stronger statements about the EOS and the behavior of the and curves.
iii.1 Current EOS constraints
The a priori constraints on the EOS that we will use are:

The EOS can be considered known below a certain density , and we use a fixed EOS below that density.

The EOS must be a monotonically increasing function (, where is the energy density) to satisfy thermodynamic stability.

The speed of sound must be less than the speed of light, ensuring causality.

The EOS must allow for a maximum NS mass greater than observed masses. Recent pulsar mass measurements in two neutron starwhite dwarf binaries have provided convincing evidence that the maximum NS mass is . The pulsar J16142230 was found to have a mass of ( confidence) Demorest et al. (2010), and the pulsar J0348+0432 was found to have a mass of ( confidence) Antoniadis et al. (2013). We will take the lower bound of on the mass of J0348+0432 as a solid lower bound on the maximum NS mass.
Although other constraints exist, such as the constraint that the maximum massshedding (Kepler) frequency be greater than observed spin frequencies, these turn out to be less useful, so we will focus on the ones listed above. See Ref. Read et al. (2009) for a discussion of other constraints.
iii.2 Choice of parameterization
We seek a parameterized EOS that will satisfy the above constraints and also have enough freedom to accurately fit the true EOS. Several choices have been presented in the literature. Read et al. examined various types of piecewise polytropes, defined below, and found that a fourparameter fit could adequately match a wide range of theoretical models Read et al. (2009). Steiner et al. used a model with four nuclear parameters around nuclear saturation density and an additional 4 parameters at higher densities to define a 2piece polytrope with variable dividing densities Steiner et al. (2010), and they also examined several variations Steiner et al. (2013). Finally, Lindblom constructed a spectral expansion of that appears to converge to tabulated EOS models with fewer parameters than the 4parameter piecewise polytrope constructed by Read et al. Lindblom (2010). All of these models can be made to satisfy the constraints in Section III.1. However, because the 4parameter piecewise polytrope Read et al. (2009) has been more commonly used in the GW literature, we will focus on it in this paper, and leave a detailed comparison of how the results depend on the choice of parameterization to future work.
For each piece of a piecewise polytrope, the pressure in the restmass density interval is defined by
(17) 
where is the adiabatic index and the constant is chosen such that is continuous at the boundary . For the core of the star, we use 3 polytropes with adiabatic indices , , and separated by the fixed dividing densities g/cm and g/cm. These fixed dividing densities were chosen to minimize the leastsquares error between the piecewise polytrope fit and a set of 34 tabulated theoretical EOSs Read et al. (2009). In addition to the three adiabatic indices, we also require an additional parameter to determine the overall pressure scaling which we choose to be , the pressure at the first dividing density. We will therefore use as our four EOS parameters .
For the lower density crust of the NS we use a 4piece polytrope fit to the SLy EOS given by Table II of Ref. Read et al. (2009)^{1}^{1}1The values for in Table II of Ref. Read et al. (2009) incorrectly give the pressure in massdensity units . The values of in the table should therefore be multiplied by ., and we note that the particular choice of the crust EOS effects the results at less than a percent level. We join the fixed crust EOS and parameterized core EOS together at the density where the two EOSs intersect. This usually occurs below , and we reject any set of EOS parameters where the joining density is below the start of the last polytrope piece for the crust at g/cm, where we consider the EOS to be known. Finally, given this EOS, the energy density , used for solving the stellar structure equations that determine the radius and tidal deformability, can be evaluated by integrating the first law of thermodynamics
(18) 
We note that for NSs around found in BNS systems, the central density is close to or just above the last dividing density g/cm for most EOS parameters Read et al. (2009). Information about therefore cannot come from BNS inspiral observations alone. However, when used in conjunction with the priors from causality and high mass NS observations, the inspiral will help constrain and the EOS above .
In Fig. 1, we show the constraints placed on the parameter space by the causality requirement and the existence of a NS with a mass of at least . Although the causality requirement must hold for all densities, we have only assumed that the piecewise polytrope used here is a sufficiently good fit to the true EOS for densities below , the central density of the maximum mass NS. Above this density the EOS could take some other form that is not well described by the expression . We have therefore used the weaker constraint of excluding EOS parameters that result in below but accepting those parameters if above as discussed in Ref. Read et al. (2009). Because most of the mass in a NS at its maximum allowed mass is above , the maximum mass and maximum speed of sound are mostly independent of , the adiabatic index below . We therefore show the constraint in the 3parameter subspace with , the value that restricts the other parameters the least.
In our analysis below, we will also use the constraints , , , and . We also use the above requirement that g/cm which restricts small values of when is large as shown in Fig. 4 of Ref. Read et al. (2009). These boundaries are large enough to incorporate the 34 EOSs considered in Ref. Read et al. (2009), and are also large enough that they have a minimal impact on the results below.
iii.3 Comparison between parameterized and theoretical EOSs
The piecewisepolytrope fit was originally constructed to match a wide range of theoretical EOSs whether or not they satisfied the then current constraints Read et al. (2009). In this paper we will examine the 7 EOSs that satisfy the causality constraint for densities up to and also have maximum masses above . We list these EOSs in Table 1 along with the piecewisepolytrope parameters that minimize the leastsquares residual Read et al. (2009). We also list several NS properties and the associated errors in reproducing them with the fit. In Fig. 2, we compare the radius and tidal deformability for these EOSs with their leastsquares fits. The error in and above for the EOS fits is usually less than 10%, except for the ALF2 EOS.
EOS  residual  %  %  %  %  

SLy  34.384  3.005  2.988  2.851  0.0020  0.989  1.41  2.049  0.02  11.736  0.21  1.69  1.10 
ENG  34.437  3.514  3.130  3.168  0.015  1.000  10.71  2.240  0.05  12.059  0.69  2.20  4.93 
MPA1  34.495  3.446  3.572  2.887  0.0081  0.994  4.91  2.461  0.16  12.473  0.26  2.78  2.47 
MS1  34.858  3.224  3.033  1.325  0.019  0.888  12.44  2.767  0.54  14.918  0.06  8.13  4.17 
MS1b  34.855  3.456  3.011  1.425  0.015  0.889  11.38  2.776  1.03  14.583  0.32  7.28  4.69 
H4  34.669  2.909  2.246  2.144  0.0028  0.685  4.52  2.032  0.85  13.774  1.34  5.12  5.40 
ALF2  34.616  4.070  2.411  1.890  0.043  0.642  1.50  2.086  5.26  13.188  3.66  4.27  24.34 
Iv Bayesian inference of EOS parameters
iv.1 Derivation
Given the GW model and EOS fit described in Sections II and III, we will describe here a method for estimating the EOS parameters from a set of GW observations. We want to find the posterior density function (PDF) for the universal EOS parameters that are common to all NSs Glendenning (1996) and the waveform parameters that are unique to each of the BNS inspiral events^{2}^{2}2The EOS parameters are related to the waveform parameters through the relations and , and it would be possible to use the EOS parameters instead of and as part of the waveform parameters. However, it is simpler to not modify existing parameter estimation codes, and treat them as separate parameters until Eq. (26).. The data is composed of the data streams from the GW detector network for each of the inspiral events. represents the waveform model, chosen here to be the TaylorF2 waveform with parameters , as well as the EOS model, chosen to be the fourparameter piecewise polytrope with parameters . represents the background information for the waveform and EOS parameters. The PDF is given by Bayes’ theorem
(19) 
The quantity is the prior probability density for the EOS and waveform parameters, is the likelihood, and the normalization constant , which we will not need to calculate here, is the evidence. We will then want to integrate over the waveform parameters to obtain a marginalized PDF for just the EOS parameters
(20) 
This dimensional integral is not easily computed, so we will decompose it into blocks and calculate it in a two step procedure. In the first step, we use an existing MCMC algorithm to sample the posterior for the parameters of each BNS event, then marginalize over the extrinsic/nuisance parameters ^{3}^{3}3 is an intrinsic parameter that in principle provides additional EOS information. However, it is unmeasurable with aLIGO Wade et al. (2014), so we treat it as a nuisance parameter and group it with the extrinsic parameters when we marginalize over them. to obtain a quasilikelihood (Eq. (25)) for the intrinsic parameters that are relevant for measuring the EOS. In the second step, we evaluate Eq. (20) by constructing a joint likelihood for the BNS events from the quasilikelihoods for each event, then reexpress in terms of the EOS parameters and masses, and marginalize over the masses in Eq. (27) using another MCMC algorithm. We obtain these expressions from Bayes’ theorem (Eq. (19)) as follows.
The prior can be decomposed into EOSparameter and waveformparameter parts using the product rule as well as the fact that the sets of waveform parameters are independent
(21) 
The conditional prior for the waveform parameters of each binary can be further decomposed with the product rule and the fact that the intrinsic and extrinsic parameters are independent
(22) 
Likewise, the likelihood for the independent BNS observations is
(23) 
and we have used the fact that the likelihood only depends on the waveform parameters to write . (The waveform signal depends on only through which is already included as a waveform parameter.)
The marginalized PDF (Eq. (20)) is now
(24) 
where we have defined the quasilikelihood for the intrinsic parameters as
(25) 
Because is a deterministic function of , and the EOS parameters,
(26) 
The marginalized PDF finally becomes
(27) 
The problem has now been reduced to computing the quasilikelihood (Eq. (25)) for each BNS event, then computing Eq. (27).
iv.2 Likelihood and signal to noise ratio
The final ingredient we need to evaluate the marginalized PDF is an expression for the likelihood for each GW event^{4}^{4}4In the following subsections, when we discuss the likelihood for individual GW events, we omit the event index for brevity.. In this paper we assume that each detector in the network has stationary, Gaussian noise and that the noise between detectors is uncorrelated. This means that the power spectral density (PSD) of the noise in detector is
(28) 
where is the Fourier transform of the noise of detector , and represents an ensemble average. For a GW event with true parameters , resulting in the GW signal , the data stream of detector will be
(29) 
For stationary, Gaussian noise, it is well known that the probability of obtaining the noise time series is
(30) 
where is the usual inner product between two time series and weighted by the PSD
(31) 
If we have a GW model that approximates the true signal , then the likelihood of obtaining the data given the GW model is
(32) 
If the model differs from the true signal , a systematic error will be introduced in the recovered waveform parameters. For a network of independent detectors described by the set of time series , the likelihood is then
(33) 
The observed signal to noise ratio (SNR) for the data from detector is a Gaussian random variable given by
(34) 
The SNR for a network of detectors is then . This is a measure of the relative power a given GW signal produces in a network of detectors.
iv.3 Averaged likelihood
In order to determine the characteristic ability of GW detectors to measure EOS parameters, we will want to average our results for Eq. (27) over many realizations of the detector noise. Ref. Nissanke et al. (2010) defined the averaged likelihood as the geometric mean of an ensemble of identical detectors measuring the same GW event
(35) 
Taking to be large, they found that this is equivalent to setting the noise used to generate the data equal to zero, but still using the characteristic noise PSD in the expression for the likelihood:
(36) 
This zeronoise likelihood will allow us to examine characteristic measurement uncertainties in the EOS parameters independent of the particular noise realization. It will also allow us in Section VI to separate the effects of systematic errors due to uncertainties in the waveform model from effects due to individual noise realizations.
iv.4 Implementation
We evaluate Eq. (27) for the marginalized PDF in a twostep procedure, similar to that described by Steiner et al. Steiner et al. (2010), where we first evaluate the quasilikelihood (Eq. (25)) for each of the BNS systems, then evaluate Eq. (27). In the first step where we evaluate the quasilikelihood , we use the MCMC sampler LALInferenceMCMC included in the LSC Algorithm Library lal as implemented in Ref. Wade et al. (2014) and described in more detail in Ref. Veitch et al. (2014). The prior for the extrinsic parameters is given in Section IIC of Ref. Veitch et al. (2014) with the additional uniform prior on of . The priors for the intrinsic parameters of interest are not contained in the quasilikelihood, so we are required to use flat priors. In most cases we used as in Ref. Veitch et al. (2014) and as in Ref. Wade et al. (2014). However, when we injected BNS systems with component masses of in Section V.3, we used and so that the posterior was minimally affected by the prior. For the likelihood , we always use the TaylorF2 waveform as our GW model because it is the fastest to generate. After running the MCMC sampler, the marginalized distribution is evaluated from the chain of samples with a Gaussian kernel density estimator.
In the second step, we sample the parameter integrand of Eq. (27) using the affineinvarient ensemble sampler emcee ForemanMackey et al. (2013). For the prior on the EOS parameters , we use uniform distributions with boundaries as described in Section III.2 and additional boundaries from the observation and causality constraints as described in Sections III.1 and III.2. For the prior on the masses , we use uniform distributions with or, when we examine systems, . The upper limit of is sufficient to include NSs for all viable EOSs. In general, the dominant cost comes from evaluating for each BNS system at each iteration, but this is sped up by precomputing the fiveparameter function on a grid and interpolating. We perform 10 runs with random initial parameters and test for convergence with the GelmanRubin diagnostic Gelman and Rubin (1992) before joining the samples. The marginalized distribution is then given by a histogram of the samples for .
Evaluating the quasilikelihood (Eq. (25)) with LALInferenceMCMC for all events in a population is very computationally expensive, so in some cases we use the Fisher matrix approximation instead. In the large SNR limit, the difference between the estimated parameters and the true parameters of the binary system obeys a Gaussian distribution Finn and Chernoff (1993). Specifically, for parameters, the likelihood is
(37) 
where is the covariance matrix, and it is given in terms of the Fisher matrix
(38) 
by the relation . In the large SNR limit, will be approximately given by the maximum likelihood. When we use the Fisher matrix approximation we will use the mass variables and instead of and and flat priors for and . The quasilikelihood (Eq. (25)) marginalized over the extrinsic parameters is simply given by the submatrix of containing the intrinsic parameters . However, when we estimate the EOS parameters in Eq. (27), we always sample the posterior with emcee.
V Results
In this Section we characterize the ability of the aLIGO–aVirgo network to measure the EOS from a population of BNS inspiral events. For the two aLIGO detectors, we use the zero detunedhigh laser power (broadband) PSD Shoemaker (2009) which represents the design sensitivity that may be achieved for the aLIGO detectors by 2019 LIGO Scientific Collaboration et al. (2013). For aVIRGO, which has higher highfrequency noise, we use the PSD fit from Ref. Manzotti and Dietz (2012) that may be achieved by 2021 LIGO Scientific Collaboration et al. (2013). We also use the TaylorF2 waveform as both the injected GW signal and the GW model used for parameter estimation.
v.1 Baseline BNS population
We start with a BNS population that has a realistic distribution of masses, number of events, and a moderate EOS, then later examine how these choices effect the results. We sample our population as follows:

Masses. NSs in most BNS systems are thought to undergo little accretion after their formation, and are therefore found to be in the relatively narrow mass range characteristic of nonaccreting NSs. The 10 currently known BNS systems Lattimer (2012) have most likely NS masses in the range – (some of these have significant uncertainties) and are shown in Fig. 2. Özel et al. Özel et al. (2012), for example, modeled the mass distribution of the known BNS systems as a Gaussian and found the most likely values for the mean and standard deviation to be and respectively. They also found that the mass ratios of the known BNS systems are consistent with each NS being drawn from this distribution independent of its companion. For simplicity, we draw the mass of each NS independent of its companion from a uniform distribution between and .

Events. Significant uncertainty exists in the BNS inspiral event rate. Ref. Abadie et al. (2010) compiled rate estimates from several populationsynthesis models and observations (Table 6 of Ref. Abadie et al. (2010)), and summarized the results as follows: a lower 95% confidence bound of 1 event per Milky Way Equivalent Galaxy per Myr (MWEG Myr), a most likely “realistic” value of 100 MWEG Myr, and an upper 95% confidence bound of 1000 MWEG Myr. This corresponds to GW detection rates of 0.4 yr, 40 yr, and 400 yr respectively for a single aLIGO interferometer, using the broadband PSD with a threshold , averaging over sky location and orientation, with NS masses of , and a density of 0.0116 MWEG Mpc Abadie et al. (2010).
We simulated a year of GW data using the “realistic” event rate above. Specifically, we calculated the event rate in a volume large enough to contain all detectable BNS events. Because inspiral events are a Poisson process, we sampled the actual number of events in a year from a Poisson distribution with this rate. We then sampled the locations of these events uniformly in the volume, the orientations uniformly on a unit sphere, and the individual NS masses uniformly in . Of these systems 120 had a network SNR 8 and 30 had a network SNR 12. We performed parameter estimation for the 20 loudest (highest SNR) sources with network SNR that ranged from 63.7 down to 13.6, integrating between the GW frequencies Hz and .

EOS. We used the piecewisepolytrope fit to the MPA1 EOS given in Table 1 which has a radius and maximum mass roughly in the middle of the range, then calculated the corresponding tidal parameter for each sampled NS. Using the fit instead of the tabulated EOS separates the systematic error due to the inexact EOS fit from the analysis of statistical errors presented here, and we leave the discussion of these systematic errors to Section VI.2.
Using zeronoise data described in Section IV.3, we show in Fig. 3 the measurability of the EOS, radius, and tidal deformability for the loudest 20 events in our population. The contours represent the 68% (), 95% (), and 99.7% () credible regions. In the left panel we plot the credible region for as well as for , where we call the piecewisepolytrope fit to the MPA1 EOS the “true” EOS. We generate these figures as follows: At the density , we evaluate for each set of EOS parameters from the MCMC simulation. We then evaluate the credible interval at that density from the sampled values. The same is done for . In the top right panel, we generated the credible interval for each mass from the radii samples that were found by evaluating for each set of EOS parameters in the MCMC simulation. For masses greater than the prior, some of the sampled EOS parameters do not allow for a stable NS. For those sampled EOS parameters, an object of that mass would lead to a black hole. In this case, the distribution of radii becomes bimodal, with a delta function at the Schwarzschild radius and weight proportional to the fraction of samples that do not allow for a stable NS. The credible interval then represents the fraction of MCMC samples that produce a NS in that radius interval or a black hole. The bottom right panel shows the confidence intervals for the tidal deformability . For the samples that produce a black hole above , Binnington and Poisson (2009).
Fig. 3 has sharp peaks in the fractional uncertainty around the variable transition density and at the fixed transition densities and between the polytrope pieces. This is due to the choice of parameterized EOS model. Allowing the transition densities to be additional free parameters would likely smooth these features out. Indeed, the EOS fit in Steiner et al. Steiner et al. (2010) which included the transition densities as free parameters did not show such features in the results for .
Above the transition density , the results increasingly underpredict the pressure. This occurs because, although the MPA1 EOS is causal with (Table 1), the corresponding piecewisepolytrope fit overpredicts by and is therefore acausal at high densities. However, the accepted MCMC samples are required to have , resulting in accepted samples corresponding to smaller pressures.
The credible interval is largest at densities below and for corresponding low mass stars where the densities are lower. This results because the bulk of NS matter is above , and we included minimal a priori information on how the core EOS joins onto the lowerdensity crust EOS. In contrast, Steiner et al. Steiner et al. (2010) parameterized the EOS around in terms of the baryon density and proton fraction with 4 free parameters (Eq. (33) of Ref. Steiner et al. (2010)), and this provides stronger a priori constraints on the behavior below (Fig. 8 of Ref. Steiner et al. (2010)). Overall, it is clear that in some density regions, a significant contribution to the credible interval comes from our choice of EOS parameterization which was not optimized for the purposes here, rather than from the sensitivity of the GW detectors.
We also find that the error in the tidal deformability is smallest in the mass interval – where the BNS masses were drawn from. This is not surprising. However, can still be measured with comparable accuracy for a much larger range of masses. This is in contrast to the results of Del Pozzo et al. Del Pozzo et al. (2013) that did not incorporate the additional information about the EOS presented in Section III.1.
Finally, we show how the credible region depends on the number of events in Fig. 4. The dashed gray curve represents the lower limit set by the priors without any BNS inspiral data. This lower limit for the radius of km above is in mild tension with Guillot et al. Guillot et al. (2013) who, combining data from observations of several NSs, found that the NS radius is 7.6–10.4km (90% confidence). A similar lower bound from the maximum mass and causality constraints was found in Ref. Hebeler et al. (2013) (dotted curve in the left panel of Fig. 11 of Ref. Hebeler et al. (2013)). However, by softening the EOS so that whenever their piecewise polytrope parameterization became acausal, they were able to weaken this constraint, and their lower limit on the radius is km less than the one presented here. The upper limit from the prior (not shown here) is a few times larger than the scale of the figures and is set mainly by the causality constraint. This is in contrast to Ref. Hebeler et al. (2013). Because they used strong assumptions about the EOS below from chiral effective field theory, they were able to place an upper limit on the NS radius of km. Because we use a much less restricted lowdensity EOS, the pressure at higher densities is allowed to have much larger values, resulting in larger radii.
When data from the BNS inspirals are included, the overwhelming majority of the information about the EOS is obtained from just the loudest 5 events in the population. After the loudest 5 events, including more events does not improve the measurability of the EOS parameters, radius, or tidal deformability.
v.2 Highest known NS mass
As discussed in Section III.1, the highest mass NSs with rigorous constraints are Demorest et al. (2010); Antoniadis et al. (2013), and we have used as the lower bound on the maximum mass. However, there is also evidence for NSs with higher masses. In particular, a class of pulsars known as blackwidows irradiate their companions and generate outflows that are accreted onto the pulsar, significantly increasing the pulsar’s mass. Using spectra to determine the radial velocity of the companion, PSR B1957+20 was found to have a mass of after correcting for the anisotropic emission of the companion which causes the center of light to lie inward relative to the center of mass. However, considering systematic uncertainties in the model they found a conservative lower limit of van Kerkwijk et al. (2011). Another blackwidow pulsar, PSR J13113430, was found to have a mass of , but considering similar uncertainties, a lower limit of was claimed Romani et al. (2012).
In Fig. 5 we demonstrate that the confirmation of a NS with a lower mass bound of , for example, would place a significantly tighter lower bound on the pressure for . The maximum mass of the MPA1 EOS is ( for the piecewisepolytrope fit), so this EOS would almost be ruled out. In contrast, the upper bound on the pressure above g/cm in Fig. 5 comes from the causality requirement. (Recall from Fig. 1, the causality constraint restricts large values of and partially restricts large values of .) Higher mass measurements will therefore not decrease the upper limit on the pressure at the highest NS densities. Similarly, observations of higher mass NSs place tighter lower bounds on the radius and tidal deformability at high masses, but do not improve the upper bound.
v.3 Distribution of BNS masses
The density profile in a NS depends on its mass, with more massive stars consisting of denser matter. We thus expect low mass stars to better estimate the lower density EOS, and higher mass stars to better estimate the higher density EOS. We generate four BNS populations with different mass distributions, then examine how the error in recovering the EOS depends on the NS masses. The first population is the same as above, using masses uniformly sampled in the range –. We also examine three additional populations where the NSs are either all , all , or all . In order to make a direct comparison between these populations, we hold all of the parameters fixed except for the masses, then adjust the tidal parameters of each system as determined by the fit to the MPA1 EOS. Additionally, we only examine the loudest 5 systems in each population which, as shown above, contain the majority of the EOS information. As a result, adjusting their masses in this range will not push these events above or below the detection threshold.
In Fig. 6, we find that when all NSs have masses of either , , or , the uncertainty in pressure is smallest around g/cm, g/cm, or g/cm respectively. We find similar results for the radius and tidal deformability; the location of the minimum uncertainty scales with the observed masses, and for the tidal deformability, the minima occur very close to the masses of the observed NSs. Interestingly, the results for masses fixed at and for a uniform distribution from – are almost identical, indicating that useful information can be found even if the range of observed masses is very small.
v.4 Comparison of LALInferenceMCMC with Fisher matrix
In addition to the results found above by evaluating the quasilikelihood (Eq. (25)) with LALInferenceMCMC, we also use the Fisher matrix approximation (Eq. (37)) to evaluate the quasilikelihood. This approach is significantly faster but only accurate in the limit of large SNR signals. For the Fisher matrix approximation, we used a single detector with the broadband PSD, but scaled the amplitude such that the SNR was equal to the network SNR for the full aLIGO–aVirgo network. The waveform was identical to the one used with LALInferenceMCMC, except we did not use the term in the nexttoleading order tidal correction containing , and thus was our only tidal parameter. Because has an insignificant impact on the waveform relative to and is unmeasurable as discussed in Ref. Wade et al. (2014), this is approximately equivalent to placing a reasonable flat prior on then marginalizing over its values as is done in the LALInferenceMCMC calculation. Furthermore, we use flat priors for and in the Fisher matrix calculation instead of the flat priors for and used in the LALInferenceMCMC calculation.
In Fig. 7, we compare the results for calculating the quasilikelihood with LALInferenceMCMC and the Fisher matrix approximation. We find that even with the differences listed above, the results are still comparable. The Fisher matrix approximation, however, slightly underestimates the uncertainty in the pressure, radius, and tidal deformability.
v.5 Dependence on EOS model
In this subsection we determine how the measurability of the EOS depends on the choice of “true” EOS. To do this, we use the same population of BNS systems as in Section V.1, then vary the EOS and corresponding tidal deformability. Here we use the tabulated EOS models listed in Table 1 instead of the leastsquares fit used above. This will allow us to examine systematic errors in Section VI.2 from the inexact EOS parameterization. For efficiency we also calculate the quasilikelihood with the Fisher matrix approximation which gives results consistent with LALInferenceMCMC as shown above.
The uncertainty in the pressure for these EOSs is shown in Fig. 8, and there are a few features to note. First, if the true EOS has a maximum speed of sound below the central density of the maximum mass NS that is close to (SLy, ENG, MPA1, MS1, MS1b), the causality constraint places a useful upper bound on the pressure estimate at densities above . However, for the H4 and ALF2 EOSs where , the causality requirement only provides a weak constraint on the highdensity EOS. Second, we note that softer EOSs (lower pressures) result in stars that are more easily compressed (smaller radii) and have higher densities. They will therefore probe higher densities. We find that for the softer EOSs SLy, ENG, and MPA1, the uncertainty in the pressure is minimized at densities , whereas for the stiffer EOSs MS1, MS1b, H4, and ALF2, the uncertainty in the pressure is minimized at densities .
The corresponding uncertainties in and for these 7 EOSs is shown in Fig. 9. The full width of the 95% credible regions is –2 km for the radius and – g cm s for the tidal deformability in the mass range –, and is roughly consistent for all EOSs. For , the credible region is smallest in the range –1.6 from which the BNS population is sampled.
v.6 Noise realizations and sampled populations
In the above results, we used the zeronoise averaged likelihood (Eq. (36)) which gives results averaged over individual realizations of the detector noise. To determine how much the recovered EOS, radius, and tidal deformability can vary with the individual noise realizations, we injected the inspiral waveforms in our population into five different sets of detector noise and recovered the parameters with LALInferenceMCMC. Fig. 10 shows the recovered EOS, radius, and tidal deformability for each of our five noise realizations. For the pressure, the lower limits are roughly the same except in the region – g/cm. This results from the fact that much of the lower bound is determined by the prior and the choice of EOS parameterization. We also show results for the zeronoise data presented in Fig. 4, and this does appear to be the average of the individual noise realizations.
Finally, we examined how much the results depend on the sampled population. We sampled five different populations from the distribution described in Section V.1. These populations had between 121 and 127 events with network SNR , and the highest network SNR event in each population was between 41 and 88. We then evaluated the credible regions for the EOS, radius, and tidal deformability using zeronoise data. The results for the loudest five events in each population are shown in Fig. 11. For a year of data with the “realistic” event rate, the reconstructed EOS as well as radius and tidal deformability are only mildly sensitive to the particular realization of the number of BNS events and source parameters.
Vi Systematic errors
vi.1 Waveform model
The presentation of statistical errors above assumed an exact waveform model and a parameterized EOS that exactly fits the true EOS. At present, however, the PN inspiral waveform is only known completely to 3.5PN order in the pointparticle terms, while the leading EOSdependent tidal term enters at the same order as 5PN pointparticle terms. Failing to include 4PN and higherorder pointparticle terms can therefore bias the recovered parameters. As demonstrated in Fisher matrix studies Favata (2014); Yagi and Yunes (2014) and in an MCMC study Wade et al. (2014), the systematic error in the tidal parameter from the current waveform uncertainty is as large as the statistical error for aLIGO.
To examine the effect that uncertainty in the waveform model has on the recovered EOS, we will use a similar analysis to these previous works. We use variations in the current PN waveforms which vary only in how the waveform phase is calculated from the energy and luminosity. The waveform variations we use, described in Section II, are the TaylorF2, TaylorT1, and TaylorT4 waveforms. For the loudest five events, we injected these three waveform variations into zeronoise data, then used the TaylorF2 waveform as the template to estimate the waveform parameters with LALInferenceMCMC. In Fig. 12 we show the bias in the recovered EOS, radius, and tidal parameter that results from the waveform uncertainty. The ordering is consistent with that of Ref. Wade et al. (2014); injecting the TaylorT1 waveform and recovering with the TaylorF2 waveform overestimates the tidal parameter while injecting the TaylorT4 waveform and recovering with the TaylorF2 waveform underestimates the tidal parameter. Likewise, injecting the TaylorT1 waveform overestimates the true pressure and radius, while the TaylorT4 waveform underestimates the true pressure and radius. Overall, failing to include the correct 4PN and higher pointparticle terms can lead to a bias that is in some cases larger than the 95% statistical uncertainty.
As in Refs. Favata (2014); Yagi and Yunes (2014); Wade et al. (2014), we focused on the uncertainties in the pointparticle description. However, other matter effects in addition to the quadrupole tidal interaction used here may also need to be accounted for. Ref. Yagi (2014), for example, calculated the correction to the PN waveform from higher multipole tidal interactions. In addition, the amplification of the tidal deformation that occurs due to resonance when the GW frequency approaches each NS’s fmode frequency also leads to a small correction Flanagan and Hinderer (2008). These effects are small, but will lead to a fractional error in the recovered parameters if not properly included.
As an alternative to the PN approximation, the effective one body (EOB) formalism, which uses various techniques to resum the PN series, may converge more rapidly to the true binary waveform. EOB waveforms have been shown to accurately reproduce numerical binary black hole (BBH) waveforms. For example, Damour et al. Damour et al. (2013) compared a recent EOB implementation with nonspinning BBH simulations of the last GW cycles and found a phase difference of radians after fitting the unknown 5PN contribution in the EOB radial potential to the numerical data. Since the tidal contribution to the waveform over this same interval is usually more than a radian, the current EOB waveform may be accurate enough for removing systematic errors due to uncertainties in the pointparticle model. Tidal interactions have also been calculated in the EOB formalism for the first few multipoles to 2PN order in the EOB radial potential Bini et al. (2012), and comparisons with numerical BNS simulations have shown that EOB waveforms are consistent with the numerical waveforms, but only after calibration of currently unknown terms Bernuzzi et al. (2012); Hotokezaka et al. (2013). Unfortunately, BNS codes are not currently as accurate as BBH codes, so waveforms calibrated with numerical BNS waveforms may still bias the recovered tidal deformability.
In our analysis we injected waveforms with zero NS spin and zero eccentricity and used a nonspinning, noneccentric waveform template to recover the parameters. While NSs in known BNS systems have dimensionless spins of , not including the spin terms in the template can lead to systematic errors in that are greater than the statistical errors if the NSs have spins of Favata (2014). Likewise, the systematic error in from not including the eccentricity terms will be greater than the statistical error if the BNS system has an initial eccentricity at 10 Hz of , and this may occur for a small fraction of BNS systems formed in dense stellar environments Favata (2014). These terms will need to be included in future studies.
Finally, in order to measure the NS EOS with BNS inspiral observations, one needs to know that one is observing a BNS inspiral. Attempting to recover the tidal parameter with a BNS waveform template, when one has actually observed a BHNS or BBH inspiral, will give erroneous results. One can positively identify a BH if its mass or spin is large enough because NSs have a more restricted range of allowed masses and spins ( to satisfy causality Rhoades and Ruffini (1974) and dimensionless spin to satisfy the Kepler constraint Cook et al. (1994)). Furthermore, NSs in known BNS systems have masses in the range shown in Fig. 2 and spins Favata (2014). Restricting their analysis to alignedspin systems, Hannam et al. Hannam et al. (2013) found that it would be difficult to distinguish between NSs and BHs in binaries except for the loudest few percent of signals. However, precession can break the massratio–spin degeneracy, and using precessing waveforms, Chatziioannou et al. Chatziioannou et al. (2014) found that one can distinguish between NSs and BHs for the majority of detected signals. Fortunately, as shown in Section V.1, the majority of EOS information comes from the loudest few events, so even if one cannot determine if the weakest signals come from a BBH, BHNS, or BNS event, this will not significantly effect one’s ability to measure the EOS.
vi.2 EOS fit
In addition to errors in the waveform model, the choice of EOS parameterization can also effect the measurement of the EOS, and if unable to sufficiently capture complex behavior in the true EOS, the parameterization will introduce systematic errors. As shown in Fig. 8, the piecewisepolytrope fit can usually reproduce the tabulated EOS to a few percent. However, at the lower and upper density regions, the fit becomes worse. In addition, for the ALF2 EOS, there is a significant change in the EOS around