Constraining the neutron star equation of state using multi-band independent measurements of radii and tidal deformabilities

TL;DR

This paper combines gravitational wave and X-ray observations using Bayesian methods to tightly constrain the neutron star equation of state, improving radius bounds over previous studies.

Contribution

It introduces a joint Bayesian analysis of multi-band measurements to derive more precise neutron star EoS constraints.

Findings

Tighter radius constraints than previous bounds.

Consistent EoS parameters with multi-band data.

Improved understanding of matter at supranuclear densities.

Abstract

Using a Bayesian approach, we combine measurements of neutron star macroscopic observables obtained by astrophysical and gravitational observations, to derive joint constraints on the equation of state (EoS) of matter at supranuclear density. In our analysis we use two sets of data: (i) the masses and tidal deformabilities measured in the binary neutron star event GW170817, detected by LIGO and Virgo; (ii) the masses and stellar radii measured from observations of nuclear bursts in accreting low-mass X-ray binaries. Using two different parametrizations of the equation of state, we compute the posteriorprobability distributions of the EoS parameters, using which we infer the posterior distribution for the radius and the mass of the two neutron stars of GW170817. The constraints we set on the radii are tighter than previous bounds.

Tables

Table 1. Table 1: Median and 90% intervals for the masses M 1 , 2 subscript 𝑀 1 2 M_{1,2} and the tidal deformabilities Λ 1 , 2 subscript Λ 1 2 \Lambda_{1,2} of the two NS observed in GW170817  Abbott  et al. ( 2018b ) and for the masses M 3 , 4 subscript 𝑀 3 4 M_{3,4} and the radii R 3 , 4 subscript 𝑅 3 4 R_{3,4} of 4U 1702-429  Nättilä  et al. ( 2017 ) and 4U 1724-307  Nättilä  et al. ( 2016 ) .

$M_1\left[M_{\odot }\right]$	${\Lambda }_1$	$M_2\left[M_{\odot }\right]$	${\Lambda }_2$
${1.46}_{-0.09}^{+0.13}$	${255}_{-171}^{+416}$	${1.26}_{-0.12}^{+0.09}$	${661}_{-375}^{+858}$
$M_3\left[M_{\odot }\right]$	$R_3$ [km]	$M_4\left[M_{\odot }\right]$	$R_4$ [km]
${1.79}_{-0.30}^{+0.20}$	${12.81}_{-0.87}^{+0.44}$	${1.48}_{-0.87}^{+0.44}$	${11.06}_{-2.06}^{+1.57}$

Equations

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

$$\Gamma(p)\simeq\exp\Biggl{\{}\sum_{k=0}^{3}\gamma_{k}\left[\log\biggl{(}\frac{p}{p_{0}}\biggr{)}\right]^{k}\Biggl{\}},$$

$$\mathbf{\theta }=({\gamma }_0,{\gamma }_1,{\gamma }_2,{\gamma }_3,p_1^c,p_2^c,p_3^c,p_4^c).$$

$$\mathcal{L}(\mathbf{d}|\mathbf{\theta })=$$

$${\mathcal{L}}^{\text{EM}}(M_3,R_3)\times {\mathcal{L}}^{\text{EM}}(M_4,R_4),$$

Full text

Introduction. The detection of the gravitational wave (GW) signal emitted in the coalescence of the binary neutron star (BNS), GW170817 Abbott et al. (2017a), offers a unique opportunity to probe the properties of matter at the extreme densities occurring in a neutron star (NS) core. This detection has stimulated a number of studies in this direction, both by the LIGO/Virgo collaboration (LVC) and by independent groups. As a result, new constraints on NS internal composition have been derived, according to which the equation of state (EoS) of nuclear matter is more likely to be soft, leading to stellar configurations with high compactness and small radii Shibata et al. (2017); Rezzolla et al. (2018); Margalit and Metzger (2017); Ruiz et al. (2018); Most et al. (2018); Annala et al. (2018); Bauswein et al. (2017); Coughlin et al. (2018a); Radice and Dai (2018); Radice et al. (2018); Abbott et al. (2018a, b); De et al. (2018); Coughlin et al. (2018b); De et al. (2018); Carson et al. (2018); Raithel et al. (2018); Zhang et al. (2018); Landry and Essick (2018); Malik et al. (2018); Tews et al. (2018); Lim and Holt (2018); Kumar and Landry (2019).

The GW signal emitted in a BNS coalescence carries the imprint of NS structure in both the inspiral and merger/post-merger phases Hinderer et al. (2010); Maselli et al. (2013); Read et al. (2013); Del Pozzo et al. (2013); Wade et al. (2014); Bauswein and Janka (2012); Bose et al. (2018). During the inspiral, the information on the stellar composition is encoded at the leading order in the quadrupolar tidal deformability $\Lambda$, which describes how the shape of one star changes in response to the external tidal field Hinderer (2008). For a given equation of state, this parameter depends solely on the stellar compactness ${\cal C}={M/R}$, i.e., on the ratio between mass and circumferential radius of the NS at equilibrium. Moreover, $\Lambda$ is a monotonic function of ${\cal C}$, and decreases as the compactness grows, i.e. for more compressible (soft) matter.

The analysis of the data of GW170817 has allowed to estimate the average parameter111$\tilde{\Lambda}$ is the actual parameter that enters at the leading order into the gravitational waveform Flanagan and Hinderer (2008); Vines and Flanagan (2013); Vines et al. (2011).

[TABLE]

where $\Lambda_{1,2}$ are the individual NS tidal deformabilities and $M_{1,2}$ are the NS masses Abbott et al. (2018a, b). Assuming low-spin priors, the analysis carried out by the LVC collaboration yields a value of $\tilde{\Lambda}=300^{+420}_{-230}$ at 90% confidence level Abbott et al. (2018a). Combined with the posterior distributions of the inferred masses, this result leads to a constraint on the NS radius of $10\lesssim R\lesssim 13$ km, which excludes EoS predicting stiff matter, i.e., less compact stars Abbott et al. (2018b).

The first BNS event also marks the dawn of multi-messenger astronomy, which will combine observations in the electromagnetic (EM) and in the gravitational bandwidths, expected to cover a broad range of wavelengths, and different stages of the evolution of the observed sources.

Following the discovery of the EM counterpart of GW170817, the information obtained from the GW data has been complemented by that inferred from the electromagnetic observations associated with the merger/post-merger phase of the binary event Abbott et al. (2017b, c); Coulter et al. (2017); in particular, the properties of the gamma ray burst and the kilonova light curves have been exploited Rezzolla et al. (2018); Margalit and Metzger (2017); Coughlin et al. (2018b); Radice and Dai (2018); Radice et al. (2018); Bauswein et al. (2017); Coughlin et al. (2018a); Most et al. (2018). Thus, in these studies EM and GW observations of the same event have been fully exploited.

In this paper we combine independent measurements of NS macroscopic observables obtained from EM and GW data to derive joint constraints on the EoS of matter at supranuclear density. We consider two distinct datasets based on: (i) masses and tidal deformabilities extracted from the the data analysis of GW170817; (ii) masses and radii measured through spectroscopic observations of NS thermonuclear bursts in low-mass X-ray binaries Nättilä et al. (2016, 2017). Using a phenomenological parametrized EoS, we infer the posterior probability distribution of the EoS parameters through a fully Bayesian analysis, and derive new bounds on the radius of the two NS coalescing in GW170817. Our results show how our understanding of the EoS of matter at supranuclear density can benefit from the synergy of data coming from astrophysical phenomena spanning very different dynamical regimes, and detected with distinct experimental setups.

Furthermore, the proposed approach is a promising tool to exploit the data of high-precision surveys which, in the near future, will be available from space satellites Arzoumanian et al. (2014) and from advanced and third generation GW interferometers Abbott et al. (2017d); Hild et al. (2010).

Parametrized EoS. The thermodynamical properties of matter inside a cold neutron star can be described by a barotropic relation between the pressure $p$ and the energy density $\epsilon$, i.e. by the equation of state, $p=p(\epsilon)$. At densities $\rho\lesssim\rho_{0}$, where $\rho_{0}\sim 2.7\times 10^{14}\mathrm{g/cm^{3}}$ is the equilibrium density of nuclear matter, the EoS has been determined by extrapolating the results of terrestrial experiments on atomic nuclei Haensel and Pichon (1994); Abrahamyan et al. (2012); Tsang et al. (2009); Tamii et al. (2011), and there is a general consensus on its properties. For densities above the saturation point, typical of a NS core, the EoS is less certain; indeed, due to the complexity of quantum chromo-dynamics and to the difficulty of testing these regimes with experiments on Earth, hadronic interactions are described by a variety of models based on different approaches and assumptions Lattimer and Prakash (2001). These EoS, when used to describe a NS, lead to different values of the observables and to different relations between radius, or tidal deformability, and mass. Thus, astrophysical and gravitational wave measurements of these observables can be exploited to constrain the EoS, solving the so-called inverse stellar problem Ozel and Psaltis (2009); Read et al. (2009a); Steiner et al. (2013); Lindblom and Indik (2012, 2014); Lackey and Wade (2015); Guillot and Rutledge (2014); Ozel et al. (2016); Raithel et al. (2017); Abdelsalhin et al. (2018); Carney et al. (2018); Lindblom (2018a). Phenomenologically parametrized EoS are particularly useful in this respect, as they allow us to describe a large class of theoretical EoS through a relatively small set of coefficients Read et al. (2009b); Lindblom (2010); Steiner et al. (2010); Raithel et al. (2016); Lindblom (2018b). Moreover, they provide a unique tool to combine stellar parameters obtained from NS observations in different waveband, to infer features of the true EoS which may not be predicted by current models.

In this work we focus on the spectral representation developed by Lindblom Lindblom (2010). This model is based on a series expansion of the adiabatic index $\Gamma(p)$. It has been shown that most theoretical EoS are well approximated including the first four terms in the expansion, which correspond to the 4 free parameters of the model $(\gamma_{0},\gamma_{1},\gamma_{2},\gamma_{3})$, such that

[TABLE]

where $p_{0}$ is the pressure at the crust-core interface, which we choose as in Abbott et al. (2018b).

The Bayesian framework. The goal of this work is to combine measurements of NS observables obtained in different astrophysical channels, such as masses, radii and tidal deformabilities, and to reconstruct the parameters of a phenomenological EoS. The mass and the radius of a NS can be measured through electromagnetic observations in the X-, optical and radio wavebands of low-mass binary systems Özel and Freire (2016); Nättilä et al. (2016, 2017). The detection of gravitational waves emitted in BNS merging allows to estimate masses and tidal deformabilities of the coalescing bodies Abbott et al. (2017a, 2018a, 2018b, 2018c). We shall now show how, using a Bayesian scheme of inference, such complementary information can be combined to put stronger constraints on the NS EoS.

Within the spectral parametrization, the equilibrium configuration of a $i=1,\ldots n$ NSs is completely specified by the $m+n$ parameters, namely by the $m=4$ coefficients $\gamma_{k}$ introduced in eq. (2), and and by $n$ values of the central pressure $p_{i=1,\ldots n}^{c}$. On the other hand, each NS observation provides $2$ observables, either the mass and the radius $(M,R)$, or the mass and the tidal deformability $(M,\Lambda)$, obtained in the EM and GW channel, respectively. Therefore, if we want to fully characterise the EoS, and infer all the spectral parameters, we need at least $N=4$ observations, yielding $2N=8$ macroscopic observables. This is the minimum number of data which would allow to determine the 8 unknown quantities

[TABLE]

Hereafter $p^{c}_{1,2}$ ($p^{c}_{3,4}$) will correspond to the pressures of 2 NSs observed in the GW (EM) channel.

The LIGO/Virgo collaboration has detected one binary NS merger so far, GW170817, and the mass and tidal deformabilities of the two stars have been estimated Abbott et al. (2018a, b). Estimates of NS masses and radii based on electromagnetic observations have been obtained using a wide variety of methodologies applied to different astrophysical environments. In this work we consider the studies carried out by Nättilä and collaborators, in which the NS masses and radii are reconstructed by using the cooling tail method for the low-mass X-ray binary 4U 1724-307 Nättilä et al. (2016) and by fitting the X-ray bursting NSs directly to the observed spectra of 4U 1702-429 Nättilä et al. (2017). GW170817 and two EM measurements provide enough information to solve the inverse stellar problem described before, with a joint set of data given by $\bm{d}=\{(M_{1},\Lambda_{1},M_{2},\Lambda_{2}),(M_{3},R_{3}),(M_{4},R_{4})\}$.

Using a Bayesian approach, we compute the posterior probability density function (PDF) of the EoS parameters given the experimental data, $\mathcal{P}(\bm{\theta}|\bm{d})\propto\mathcal{L}(\bm{d}|\bm{\theta})\mathcal{P}_{0}(\bm{\theta})$, where $\mathcal{P}_{0}(\bm{\theta})$ is the prior on the parameters, and $\mathcal{L}(\bm{d}|\bm{\theta})$ is the likelihood function. The latter, in our case, reads:

[TABLE]

where $\mathcal{L}^{GW}(M_{1},\Lambda_{1},M_{2},\Lambda_{2})$ and $\mathcal{L}^{EM}(M_{i},R_{i})_{i=3,4}$ are the probability computed by the LIGO/Virgo collaboration Com ; Ind and and Nättilä et al. Nättilä et al. (2016, 2017), respectively. We sample the posterior distribution using Markov chain Monte Carlo (MCMC) simulations based on the Metropolis-Hastings algorithm Gilks et al. (1996). The MCMC convergence is enhanced by a Gaussian adaptation algorithm Kjellstrom and Taxen (1981); Müller and Sbalzarini (2010) (see Abdelsalhin et al. (2018) and reference therein for a detailed discussion on this approach).

Numerical setup. For the GW data, the likelihood $\mathcal{L}^{GW}(M_{1},\Lambda_{1},M_{2},\Lambda_{2})$ is given by the joint probability distribution for the masses and tidal deformabilities of GW170817 inferred by the LVC using a parametrised EoS Abbott et al. (2018b). For the EM sector, $\mathcal{L}^{EM}(M_{i},R_{i})_{i=3,4}$ is given by the joint probability distributions of the most accurate measurements of masses and radii provided in Nättilä et al. (2016, 2017). The latter correspond to two NSs in low-mass X-ray binaries observed during thermonuclear bursts, namely 4U 1724-307 Nättilä et al. (2016) and 4U 1702-429 Nättilä et al. (2017). Table 1 shows the median and the 90% confidence intervals for the GW and EM data.

We assume a flat prior distributions for all parameters, namely: $\gamma_{0}\in[0.2,2]$, $\gamma_{1}\in[-1.6,1.7]$, $\gamma_{2}\in[-0.6,0.6]$ and $\gamma_{3}\in[-0.02,0.02]$, in agreement with the analysis made in Ref. Abbott et al. (2018b). Furthermore, for each combination of parameters, we require the corresponding adiabatic index to be $\Gamma(p)\leq 7$ up to the central pressure of the NS maximum mass configuration. Finally, we ask all central pressures to be uniformly distributed in the range $p^{c}_{i}\in[10^{33},10^{38}]$ dyn/cm2. Also, we ask that the inferred EoS is consistent with existent astrophysical and theoretical bounds. Specifically, we require that: (i) each EoS model supports a maximum mass $M_{max}\geq 1.97M_{\odot}$ Antoniadis et al. (2013), and (ii) the EoS remains causal (speed of sound $c_{s}\leq c$, where $c$ is the speed of light in vacuum) up to the maximum mass configuration. For both models, we run 24 independent chains of $n=2\cdot 10^{7}$ samples, discarding the first $\sim 10\%$ points as burn-in. The convergence of the MCMC is determined by monitoring each chain autocorrelation function and cross-checking the chains through a standard Rubin test Gilks et al. (1996).

Results. The four panels of Fig 1 show the posterior distributions for the parameters of the spectral equation of state. In each panel, dashed vertical lines identify 90% confidence intervals.

The lowest order term in the expansion of the adiabatic index $\Gamma(p)$, i.e. $\gamma_{0}$, provides the weakest bounds with respect to the priors, while the linear, $\gamma_{1}$, and the “high-density”, $\gamma_{2}$ and $\gamma_{3}$, coefficients lead to the strongest constraints. Such posteriors are consistent with previous studies on large density regimes within NS cores featuring quark matter Most et al. (2018); Annala et al. (2018). As indicated in (3), the central pressures are also sampled in the MCMC. We find that for the neutron stars of GW170817 these quantities are in agreement with the values quoted by the LVC Abbott et al. (2018b).

Having inferred the probability distributions for the parameters (3), we can compute the posteriors for radius and mass of the two NS of GW170817, by solving the TOV equations for each of the model sampled by the MCMC.

Figure 2 shows the joint distribution of mass and radius of the two stars of GW170817, together with the confidence intervals derived by LVC Abbott et al. (2018b). The $M-R$ posteriors derived with our approach fit within the LIGO/Virgo distributions inferred through GW data alone. The values of marginalised mass and radius for the two stars, at 90% confidence level are $M_{1}=1.45^{+0.08}_{-0.06}\ M_{\odot}$, $R_{1}=12.36^{+0.52}_{-0.38}$ km and $M_{2}=1.28^{+0.05}_{-0.06}\ M_{\odot}$, $R_{2}=12.32^{+0.66}_{-0.43}$ km. While the intervals of reconstructed masses are close to the values inferred by LVC, the posterior distributions along the $R-$direction obtained by the multi-messenger approach are now effectively narrower.

This is more evident in Fig. 3 in which we plot, for the same GW170817 stars, the corresponding marginalised distributions of the two radii. The dashed and the solid vertical lines identify the 90% credible interval determined by LVC Abbott et al. (2018b), and by our approach, respectively. As already seen in Fig. 1, the analysis we perform, which combines GW and EM observations, is consistent with the results derived by LIGO/Virgo alone. Most notably, the inclusion of the new datasets shrinks the posterior distribution of $R$, which is now determined to an accuracy below $10\%$. Note that including in our analysis a more massive star, i.e. the NS with mass and radius estimated through EM observations (specifically $M_{3}$ in Table 1), we are probing the EoS in a region where the energy density is larger with respect to that probed by the LVC analysis.

Conclusions. Multi-wavelength observations of relativistic sources provide an arena where the joint efforts of the astrophysics, high-energy and particle physics community convey to provide new insights on the fundamental laws of Nature. Neutron stars are among the primary targets of this quest, as unique laboratories to investigate the behaviour of matter at densities not reproducible in experiments on Earth. The detection of the first coalescing binary composed of two neutron stars has allowed the LIGO/Virgo collaboration to derive the first GW constraint on the equation of state of matter in the inner core of neutron stars. A large variety of follow-up analyses have been pursued to further exploit the observation of the electromagnetic counterpart in coincidence with the gravitational event. As a result, new bounds have been derived, which indicate that the EoS of matter in the inner core of a neutron star is in the soft sector, and produces more compact stars.

In this paper we have made a step forward in this search, by combining independent measurements of NS macroscopic parameters, namely the radius and the tidal deformabilities from the LIGO/Virgo event, and the mass and the radius derived from EM observations of low mass X-ray binaries, in the spirit of multi-messenger astrophysics. Using this approach, we have been able to set tighter constraints on the radius of the two neutron stars coalescing in GW170817, thus supporting, and strengthening, the observational evidence that neutron star cores are composed of soft nuclear matter.

Further detections by ground-based interferometers with higher sensitivities will provide new and more accurate observations, which can be used to improve the constraints presented in this Letter. It is worth remarking that phenomenological parametrizations can introduce systematics which may affect the final results of the bayesian analysis. In this regard, a wider sample of NS will also allow to test the nature and the relevance of such systematics introduced by the specific parametrization Fasano et al. .