The composition and structure of neutron stars (NSs) 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 [1,7,[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 Λ 1 and Λ 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 Λ. 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 NS with a SLy [2] equation of state, M crust /M ≈ 1% below ρ crust ≈ 1.4 × 10 14 g cm -3 ), but contribute to a larger portion of the radius ( R crust /R ≈ 6% ). For this reason, since different crust choices are also not anticipated to strongly affect the determination of the tidal deformability parameters Λ, it is not expected to be able to distinguish between different low-density models through GW analysis alone. This, in turn, may de-facto induce a 'systematic error' on the obtained measurements of the NSs radii. 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 ρ ≈ 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 Λ and M distributions, but can shift the implied radius up by ≈ 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.
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.
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 ≈ 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,43,45,46].
In general, whatever the combination of models chosen to describe the complete low-density section is, it 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 δ = 1 -2x, through a Taylor expansion around x = 1/2 of the energy density per nucleon E one finds:
with
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 + . . . ).
From [41] we retrieve a set of crust EOSs computed through a CLDM + BPS model, which reach up to approximately ρ 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 [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 1). 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.
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. 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. Therefore, we choose employ this parameterization in our PE runs.
The basic idea behind the spectral decomposition of [44] is that of expressing the adiabatic index Γ( 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 functionsf i = [ln( p/p 0 )] i , i ∈ N in our specific case-multiplied by some coefficients γ i :
Obtaining the expression of the energy density as a function of pressure requires then a simple integration:
in which e 0 , p 0 is the starting point of the decomposition in the energy density-pressure plane and
Running PE on the real data of GW170817 using spectral decomposition means that every spectral coefficient γ i , i ∈ [0, 3] is sampled, in place of the tidal deformabilities Λ 1 and Λ 2 , by stochastically walking through the parameter space [40]. In order to cover a wide range of candidate EOSs, we sample γ 0 ∈ [0.2, 2] , γ 1 ∈ [-1.6, 1.7], γ 2 ∈ [-0.6, 0.6] and γ 3 ∈ [-0.02, 0.02] . We additionally impose that the adiabatic index Γ( p) ∈ [0.6, 4.5] , as done in [1]. Each set {γ 0 , . . . , γ 3 } i can be mapped through (4) in an EOS p i(e), i = 1, . . . , N samples . The pressure-density credible levels can be obtained by choosing a list of N energy density values e j , j = 1, . . . , N, creating a pressure histogram for each value on the list by evaluating p i (ej) ∀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 Λ, 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 1), and the crust-core transition point e t, pt changed accordingly. The approximant used, IMRPhenomPv2NRTidal [54][55][56][57] (see the references contained in [1] for more detail), 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.3 -1.4 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 Λ that we would get with our new crust. We replace the SLy EOS with the (S 0 = 30, L = 65) MeV or the (S 0 = 34, L = 35) MeV crust and glue the final points of the new EOS to the final point of the original SLy crust table ( e m, pm). Then, for pressures p and energy densities e higher than p m and e m we compute e( p) through (4), using as {γ 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 Λ (figures 3(a) and (b)). The radii values found are R (30,65) pr
km and R (30,65) pr 2 = 11.7 +1.4 -1.3 km; R (34,35) pr 1 = 12.0 +1.5 -1.4 km and R (34,35) pr 2 = 12.0 +1.5 -1.4 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 ( ρ 0 = 2.8 × 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 ρ 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 [58]. Such correlations would come through S 0 and L choices rather than from the crust densities themselves.
In figure 3(a) we plot the marginalized two-dimensional posterior distributions of Λ 1 versus Λ 2 , 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 versus R marginalized two-dimensional 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 ∆R -= R m SLy -R m (30,65) and ∆R + = R m (34,35) -R m SLy . We find ∆R + 1 = 0.1 km, ∆R -1 = 0.2 km and ∆R + 2 = 0.1 km, ∆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 Λ(R) curves (figures 4(a) and (b)). Λ(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 ∆R . After defining R X Y (M) as the R(M) relation obtained by inverting the X = 5th or X = 95th M(R) percentile curve computed using the Y crust, we have that ∆R = ∆R 0 + ∆R + + ∆R -, where ∆R 0 = R 95 SLy -R 5 SLy is the original uncertainty, obtained when considering the SLy crust only, and ∆R + = R 95 (34,35) -R 95 SLy , ∆R -= R 5 SLy -R 5 (30,65) are the corrections which account for the different crusts. For M 1 = 1.57M and M 2 = 1.20M we find ∆R + 1 = 0.2 km, ∆R -1 = 0.1 km and ∆R + 2 = 0.2 km, ∆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.
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 confirm our expectation of 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. For each percentile, 5th or 95th, the Λ(M) curves are on top of each other, while the M(R) curves differ more as the mass decreases. This can be easily explained considering that for low-mass stars the crust has more relative importance than for heavier NSs. Note also how our prediction (orange dashed) gives a good approximation of the curve obtained with the actual run (green).