The impact of the crust equation of state on the analysis of GW170817

TL;DR

This paper investigates how different crust equations of state affect neutron star radius estimates from GW170817 data, finding that crust choice influences radius but not mass or tidal deformability, with a systematic error of about 0.3 km.

Contribution

It demonstrates that crust EOS variations have a limited impact on mass and deformability but significantly affect radius estimates in gravitational wave analyses.

Findings

Crust EOS choice does not significantly alter mass or tidal deformability.

Different crust models cannot be distinguished through GW data.

Crust-related systematic error on radius is approximately 0.3 km.

Abstract

The detection of GW170817, the first neutron star-neutron star merger observed by Advanced LIGO and Virgo, and its following analyses represent the first contributions of gravitational wave (GW) data to understanding dense matter. Parameterizing the high density section of the equation of state (EOS) of both neutron stars through spectral decomposition, and imposing a lower limit on the maximum mass value, led to an estimate of the stars' radii of $R_1 = 11.9_{- 1.4}^{+ 1.4}$ km and $R_2 = 11.9_{- 1.4}^{+ 1.4}$ km. These values do not, however, take into account any uncertainty owed to the choice of the crust low-density EOS, which was fixed to reproduce the SLy EOS model. We here re-analyze GW170817 data and establish that different crust models do not strongly impact the mass or tidal deformability of a neutron star: it is impossible to distinguish between low-density models with GW analysis. However, the crust does have an effect on the inferred radius. We predict the systematic error due to this effect using neutron star structure equations, and compare the prediction to results from full parameter estimation runs. For GW170817, this systematic error affects the radius estimate by 0.3 km, approximately $3\%$ of the NS radii.

Figures

Tables

Table 1. Table 1: Left: Table containing the radii values obtained from the 2D M vs R distributions (figure 3 (b)). Right: Table containing the R i X subscript superscript 𝑅 𝑋 𝑖 R^{X}_{i} radii values, obtained by matching different crusts to the X = 5 t ​ h 𝑋 superscript 5 𝑡 ℎ X=5^{th} or X = 95 t ​ h 𝑋 superscript 95 𝑡 ℎ X=95^{th} pressure-density percentile curves, inverting the implied M(R) relation and fixing M 𝑀 M to M 1 = 1.57 ​ M ⊙ subscript 𝑀 1 1.57 subscript 𝑀 direct-product M_{1}=1.57M_{\odot} or M 2 = 1.20 ​ M ⊙ subscript 𝑀 2 1.20 subscript 𝑀 direct-product M_{2}=1.20M_{\odot} (figure 4 (b)).

	$R_1\left[km\right]$	$R_2\left[km\right]$	$R_1^5\left[km\right]$	$R_1^{95}\left[km\right]$	$R_2^5\left[km\right]$	$R_2^{95}\left[km\right]$
SLy (32-46) run	${11.9}_{-1.4}^{+1.4}$	${11.9}_{-1.4}^{+1.4}$	$10.5$	$13.5$	$10.5$	$13.5$
30-65 run	${11.7}_{-1.4}^{+1.4}$	${11.7}_{-1.3}^{+1.4}$	$10.4$	$13.3$	$10.4$	$13.2$
30-65 prediction	${11.7}_{-1.3}^{+1.4}$	${11.7}_{-1.3}^{+1.4}$	$10.4$	$13.3$	$10.4$	$13.2$
34-35 prediction	${12.0}_{-1.4}^{+1.5}$	${12.0}_{-1.4}^{+1.5}$	$10.5$	$13.7$	$10.7$	$13.7$

Equations

$$E(n,x)=E_0(n,1/2)+S(n,1/2){\delta }^2+\dots$$

$$S(n,1/2)=S_0+L(n-n_0)/3n_0+\dots$$

$$\Gamma (p)=exp[\sum _{i=0}^{\infty }{\gamma }_i{f}^i(p)]$$

$$e(p)=\frac{e_0}{\mu (p)}+\frac{1}{\mu (p)}{\int }_{p_0}^p\frac{\mu (p^{\prime })}{\Gamma (p^{\prime })}d{p}^{\prime }$$

$$\mu (p)=exp[-{\int }_{p_0}^p\frac{d{p}^{\prime }}{p^{\prime }\Gamma (p^{\prime })}]$$

Full text

1 Introduction

The composition and structure of neutron stars (NS) is a longstanding question for discussion in the scientific community. Knowing the properties of these very dense objects would contribute to the understanding of matter under extreme conditions, with implications for both astrophysics and nuclear physics, as reviewed in e.g. [3, 4, 5, 6].

In August 2017, the first observation of a NS-NS merger, GW170817 [7], was made by the Advanced LIGO [8] and Virgo [9] instruments. Its subsequent analysis represents the first contribution of gravitational-wave (GW) data to understanding dense matter and the NS radius $R$ [7, 1, 10, 11, 12, 13, 14, 15]. Rapid detection and localization of the signal [16, 17, 18, 19] enabled multimessenger followups [20, 21, 22, 23, 24, 25, 26], which themselves suggest various scenarios that may also limit the properties of cold dense matter, e.g. [27, 28, 29, 30]. We here focus on what GW data alone implies.

Comparing GW data to solutions of the relativistic two-body problem provides a measurement of the two stars’ masses $M_{1}$ and $M_{2}$ and dimensionless tidal deformabilities $\Lambda_{1}$ and $\Lambda_{2}$, which describe the ratio of the body’s tidally induced quadrupolar deformation to the tidal potential caused by its companion, through the mass-weighted sum $\tilde{\Lambda}$. The tidal parameters depend on the compactness of the star $C=R/M$, both explicitly and through the relativistic tidal Love number $k_{2}$ [31, 32, 33, 34, 35, 36, 37, 38].

When determining $R$ from GW analysis, an added level of uncertainty comes from the choice of the crust structure model. The outer low-density layers of the star contain a small fraction of the mass (for a M = 1.4 $M_{\odot}$ NS with a SLy equation of state, $M_{crust}/M\approx 1\%$ below $\rho_{crust}\approx 1.4\times 10^{14}g/cm^{3}$), but contributes a larger portion of the radius ($R_{crust}/R\approx 6\%$). If the tidal deformability parameters $\Lambda$ are mostly unaffected by the choice of the crust, then it will not be possible to distinguish between different low-density models through GW analysis, de-facto adding a “systematic error” on the values of the NS radii thus obtained. Notably, one implication of the low sensitivity of tidal parameters to the crust densities would be that GW measurements give more direct information on higher densities, and that therefore in this region the constraints obtained from analyses are independent of uncertainties in crust.

The aim of this paper is to quantify the effect of the choice of the crust equation of state (EOS) on parameter estimation (PE) for GW170817, especially on the radius estimate of [1]. In order to do so, we re-analyse the GW data with LALInference [39], following the parameterized EOS method of [40], but modifying the low-density region. In the LIGO-Virgo analysis [1], densities below $\rho\approx 10^{14}g/cm^{3}$ were fixed to the SLy description of [2]. We replace the fixed region with crust EOSs described in [41], which were obtained through the combination of Compressible Liquid Drop Model (CLDM) [42] and Baym, Pethick and Sunderland (BPS) models [43], while continuing to parameterize the higher-density core with a spectral decomposition following [44]. In parallel, we try to predict the effect of this change without full reanalysis: by gluing different crusts to core EOSs recovered in [1], we obtain a quantitative prediction of the impact of the crust on the radii recovered for the NSs involved in the GW170817 coalescence. We find that varying the crust has negligible impact on the $\Lambda$ and M distributions, but can shift the implied radius up by $\approx 0.3$ km in the full PE.

The paper is organized as follows: section 2 is dedicated to the NSs’ EOS, focusing on its parametrization in the core of the star and on choosing appropriate crust models with respect to the existing bounds on symmetry energy and its slope, two important nuclear parameters whose meaning will be briefly presented; section 3 describes the methods used to estimate and predict the stellar parameters’ values; finally, in sections 4 and 5 the estimate of the systematic error entailed by ignoring crust variations is computed and discussed.

2 Equation of State

Statements on the behaviour of NS matter quantitatively translate into imposing constraints on its energy density - pressure relationship, the EOS, which then lead to limits on stellar parameters such as maximum mass and radius. In this section we first give a quick overview of the crust composition and of the model EOS chosen to describe it. We then introduce the spectral decomposition parametrization of the high-density EOS adopted in our LALInference run.

2.1 Crust Equation of State

NSs are objects so dense that it is possible for them to develop a sturdy crust made of very neutron-rich nuclei even at the incredibly high temperatures of their surface ($T\approx 10^{5}-10^{6}K$)[45]. Qualitatively, one could imagine that when moving farther away from the core of the star - as the temperature and density decrease - the Coulomb interactions between particles become more and more important with respect to their thermal and quantum energy. At one point, this translates into the formation of neutron-rich nuclei, their locking from nucleon plasma into a lattice and the creation of a solid layer. In the innermost part of this solid crust, still very close to the core, the extreme density conditions may cause the lattice nuclei to change shape: no longer spherical, the minimum-energy structures could vaguely resemble pasta forms. At lower densities the nuclei go back to their spherical form, but are still so neutron-rich that some neutrons “drip out” of them and, if the temperature is below a critical value, form a superfluid neutron vapour. Approaching the exterior of the crust, the drip phenomenon stops; this transition divides the inner from the outer crust, which is characterized by the presence of a heavy nuclei lattice immersed in an electron gas. To quantitatively describe this complex behaviour, a number of models have been developed for both the inner (Thomas-Fermi, CLDM) and outer crust (BPS) [42, 46, 45, 43].

In general, whatever combination of models is chosen to describe the complete low-density section can be characterized by two parameters defined at nuclear saturation density $n_{0}$: the symmetry energy $S_{0}$, which encodes the energy cost of making NS matter more neutron-rich, and its slope $L$. Following common notations and defining the fraction of neutrons to all baryonic matter as $x$ and the isospin asymmetry as $\delta=1-2x$, through a Taylor expansion around $x=1/2$ of the energy density per nucleon $E$ one finds:

[TABLE]

with

[TABLE]

The slope parameter $L$ plays an exceptionally important role in NS structure, as it is closely related to the pressure of purely neutron matter at sub-saturation densities through $p(n,1)=n^{2}/(3n_{0})(L+\dots)$.

From [41] we retrieve a set of crust EOSs computed through a CLDM+BPS model, which reach up to approximately $\rho_{crust}=10^{14}g/cm^{3}$ and cover a wide region of the $S_{0}-L$ plane. Both $S_{0}$ and $L$ have been studied and constrained by a number of independent terrestrial experiments, which performed measurements of giant dipole resonances and dipole polarizabilities, nuclear masses, flows in heavy-ion collisions and neutron-skin thicknesses [47, 48, 49]. Taking into account all constraints, $S_{0}$ should range from 30 to 32 MeV and $L$ from 40 to 60 MeV for $L$ [50]. Since, however, the acceptable ranges of $S_{0}$ and $L$ are still uncertain, we focus on the larger intervals 30 to 34 MeV for $S_{0}$ and 30 to 70 MeV for $L$ [51]. We then select the upper and lower limits EOS curves, whose parameters are respectively $S_{0}=34$ MeV, $L=35$ MeV and $S_{0}=30$ MeV, $L=65$ MeV. These should give the largest impact on neutron star structure (see figure 2). We note that, although we use variation of $S_{0}$ and $L$ to establish a realistic range of crusts, we do not enforce an extrapolation to higher densities that is consistent with the chosen parameters, but allow the core EOS to vary independently of the chosen crust.

2.2 Parametrization model

There are two necessary conditions that a NS EOS has to satisfy to be physically consistent: sound has to propagate through the star slower than light in vacuum (causality) and pressure $p$ must be a monotonically increasing function of energy density $e$. When parametrizing an EOS, it would then be particularly convenient to choose a model that automatically fulfils at least the second condition. Such a parametrisation could be obtained through piece-wise polytropes [52, 53] or by applying a spectral decomposition [44] on a basis of differentiable functions. Additionally, when fitting a known EOS, spectral fits frequently have smaller residuals than piecewise-polytrope fits, even when they are performed employing fewer parameters than piecewise-polytrope fits.

The basic idea behind the spectral decomposition of [44] is that of expressing the adiabatic index $\Gamma(p)=[(e+p)/p]dp/de$, unique up to an integration constant for each EOS, as the exponential of the sum of some smooth basis functions - $f^{i}=[ln(p/p_{0})]^{i},i\in\mathbf{N}$ in our specific case - multiplied by some coefficients $\gamma_{i}$:

[TABLE]

Obtaining the expression of the energy density as a function of pressure requires then a simple integration:

[TABLE]

in which $e_{0}$, $p_{0}$ is the starting point of the decomposition in the energy density-pressure plane and

[TABLE]

3 Parameter Estimation

Running PE on the real data of GW170817 using spectral decomposition means that every spectral coefficient $\gamma_{i},i\in[0,3]$ is sampled, in place of the tidal deformabilities $\Lambda_{1}$ and $\Lambda_{2}$, by stochastically walking through the parameter space [40]. In order to cover a wide range of candidate EOSs, we sample $\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]$. We additionally impose that the adiabatic index $\Gamma(p)\in[0.6,4.5]$, as done in [1]. Each set $\{\gamma_{0},\dots,\gamma_{3}\}_{i}$ can be mapped through (4) in an EOS $p_{i}(e)$, $i=1,\dots,N_{samples}$. The pressure-density credible levels can be obtained by choosing a list of $N^{\prime}$ energy density values ${e_{j}}$, $j=1,\dots,N^{\prime}$, creating a pressure histogram for each value on the list by evaluating $p_{i}(e_{j})$ $\forall i$, and finally finding the chosen percentiles from every histogram. To then go from an EOS of the form $e(p)$ to the determination of the stellar parameters $M$, $R$ and $\Lambda$, one has to integrate TOV equations and solve the inner-outer matching problem related to the relativistic Love number $k_{2}$. We used the publicly available code of LALSimulation[39] for such operations. By additionally imposing a lower limit of 1.97 solar masses on the maximum mass value supported by the EOS, and fixing the low-density section of the EOS to reproduce the SLy model, the radii of the two NSs involved in the binary coalescence of GW170817 were estimated to be $R_{1}=11.9^{+1.4}_{-1.4}$ km and $R_{2}=11.9^{+1.4}_{-1.4}$ km. In our LALInference Markov chain Monte Carlo (MCMC) run, the hard-coded SLy crust was switched to the $S_{0}=30$ MeV $L=65$ MeV EOS (figure 2), and the crust-core transition point $e_{t},p_{t}$ changed accordingly. The approximant used, IMRPhenomPNRTidal [54, 55, 56, 57, 58], and the choice of the other priors match the settings outlined in [1]. The radii found then are $R_{1}=11.7_{-1.4}^{+1.4}$ km, $R_{2}=11.7_{-1.4}^{+1.3}$ km. Pressure-density credible levels have been computed, and are shown in figure 2.

While full MCMC runs can paint a precise picture of what happens when changing the low-density EOS, they also are very computationally-expensive and time demanding. For this reason, it is worth trying to make some rough - but fast - predictions, which would also allow one to check on the possible impact of some future variations, with a better idea of how well it will reflect the full analysis. Indeed, working under the assumption that both the masses and the core EOS are weakly affected by the choice of the crust, we can use the mass and spectral coefficients posteriors from [1] to predict the posterior distributions of the stellar parameters $R$ and $\Lambda$ that we would get with our new crust. We replace the SLy EOS with the 30-65 (34-35) crust and glue the final points of the new EOS to the final point of original the SLy crust table ($e_{m},p_{m}$). Then, for pressures $p$ and energy densities $e$ higher than $p_{m}$ and $e_{m}$ we compute $e(p)$ through (4), using as $\{\gamma_{i}\}$ those obtained from the posterior distribution of the previous analysis. Coupling these relations to the mass posteriors, we can compute the predicted distributions of the radii $R^{pr}$ of the NSs and of their tidal deformabilities $\Lambda$ (figures 3 (a) and (b)). The radii values found are $R_{1}^{30-65pr}=11.7_{-1.3}^{+1.4}$ km and $R_{2}^{30-65pr}=11.7_{-1.3}^{+1.4}$ km; $R_{1}^{34-35pr}=12.0_{-1.4}^{+1.5}$ km and $R_{2}^{34-35pr}=12.0_{-1.4}^{+1.5}$ km.

We again note that at densities from approximately $10^{14}g/cm^{3}$ and up the EOS curves obtained through spectral decomposition do not necessarily have the characteristics at saturation density ($\rho_{0}=2.8\times 10^{14}g/cm^{3}$) implied by the lower-density crust. Our aim, instead, is to compare directly with the results of [1], and estimate the effect that changing only the previous hard-coded outer crust has on PE. If we were to impose consistency on $S_{0}$ and $L$ when sampling the spectral coefficients, or if we extended the fixed EOS region from crust through to $\rho_{0}$ with the same $S_{0}$ and $L$, we expect that we would find increased correlation between crust and radius results as was seen in [59]. Such correlations would come through $S_{0}$ and $L$ choices rather than from the crust densities themselves.

4 Systematic Error Estimate

In figure 3 (a) we plot the 2D $\Lambda_{1}$ vs $\Lambda_{2}$ distributions, obtained through LALInference runs and through the predictions we made. They are all almost perfectly superimposed: this confirms that we are not able to distinguish between different low-density models through GW analysis alone. The choice of the crust does, however, impact the radii of the neutron stars involved in the coalescence. To get an estimate of the variation of the radii due to the choice of the crust alone, we compute the $M$ vs $R$ 2D distributions from the posteriors, as described in section 3 (figure 3 (b)). Our prediction of the median of the radii distribution $R^{m}$ for the 30-65 crust is slightly shifted with respect to the median obtained with the original SLy crust, and in excellent agreement with the actual MCMC result (table 1). We then estimate the systematic error due to the crust as $\Delta R^{-}=R_{SLy}^{m}-R_{30-65}^{m}$ and $\Delta R^{+}=R_{34-35}^{m}-R_{SLy}^{m}$. We find $\Delta R_{1}^{+}=0.1$ km, $\Delta R_{1}^{-}=0.2$ km and $\Delta R_{2}^{+}=0.1$ km, $\Delta R_{2}^{-}=0.2$ km. The radii of the NSs of GW170817 then become $R_{1}=(11.9^{+1.4+0.1}_{-1.4-0.2})$ km and $R_{2}=(11.9^{+1.4+0.1}_{-1.4-0.2})$ km.

More physical insight can be obtained by mapping the $90\%$ pressure-density credible levels curves, appropriately glued to the selected crusts, into M(R) and $\Lambda(R)$ curves (figures 4 (a) and (b)). $\Lambda(M)$ curves are indistinguishable, as expected, and once a mass value $M$ has been fixed we can retrieve from the M(R) relations an estimate of the uncertainty on radius $\Delta R$. Defining $R^{X}_{i}$ as the radius value obtained by inverting the $X=5^{th}$ or $X=95^{th}$ M(R) percentile curve and setting $M=M_{1}$ or $M=M_{2}$, we have that $\Delta R$ = $\Delta R_{0}+\Delta R^{+}+\Delta R^{-}$, where $\Delta R_{0}$ = $R^{95}_{SLy}-R^{5}_{SLy}$ is the original uncertainty, obtained when considering the SLy crust only, and $\Delta R^{+}=R^{95}_{34-35}-R^{95}_{SLy}$, $\Delta R^{-}=R^{5}_{SLy}-R^{5}_{30-65}$ are the corrections which account for the different crusts. For $M_{1}=1.57M_{\odot}$ and $M_{2}=1.20M_{\odot}$ we find $\Delta R_{1}^{+}=0.2$ km, $\Delta R_{1}^{-}=0.1$ km and $\Delta R_{2}^{+}=0.2$ km, $\Delta R_{2}^{-}=0.1$ km. This estimate, while not very different from the one obtained earlier, gives a less complete picture of the situation as pressure-density credible levels (CLs) do not map directly into mass-radius CLs. Nonetheless, it does a reasonable job in the mass range of the binary components, and suggests that the corrections become larger as the mass becomes smaller. This behaviour can be easily explained: the less massive the star, the higher the contribution of the crust to the total mass, and the bigger the radii differences owed exclusively to the arbitrary choice of the outer layers.

5 Conclusions

After considering a range of realistic crust EOSs, whose parameters $S_{0}$ and $L$ at saturation density span the intervals $[30-34]$ MeV and $[30-70]$ MeV respectively, we selected the $S_{0}=30$ MeV and $L=65$ MeV crust model and re-analyzed the data of GW170817. In parallel, we successfully predicted the outcome of the re-analysis, i.e the NS radii values. Such values were then used to compute the systematic error to be added to the radii estimates of GW170817, which amounts to a total of 0.3 km, approximately $3\%$ of R. The simple methods implemented to make predictions will likely be useful to quickly quantify the impact of the crust EOS on future radius estimates obtained through GW analyses. Finally, we find low sensitivity of tidal parameters to the EOS at lower crust densities. GW measurements give direct information on the high density EOS, independent of uncertainties in the crust.

Acknowledgments

The authors thank: John Douglas Veitch, Carl-Johan Haster and Sebastian Khan for their assistance with LALSuite installation and debugging; Benjamin Lackey and Katerina Chatziioannou for writing and sharing the plotting scripts we built on; Wynn Ho and Wolfgang Tichy for precious discussions; William Newton for supplying the low-density equations of state. The authors also thank INFN and NSF, who supported this work through the INFN-NSF/LIGO summer student exchange program, Giancarlo Trinchero, for support at the beginning of the exchange, and Alessandro Nagar, for discussions and comments on the manuscript. Finally, the authors are grateful for computational resources provided by the LIGO Laboratory and supported by the National Science Foundation Grants PHY-0757058 and PHY-0823459.