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

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.


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].
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 equation of state, M crust /M ≈ 1% below ρ crust ≈ 1.4 × 10 14 g/cm 3 ), but contributes a larger portion of the radius (R crust /R ≈ 6%). If the tidal deformability parameters Λ 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 ρ ≈ 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.

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.

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 ≈ 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 neutronrich 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 combination of models is chosen to describe the complete lowdensity 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 δ = 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 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 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.

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 Γ(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 ∈ 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

Parameter Estimation
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. 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 , 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 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 30-65 (34)(35)  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 [59]. Such correlations would come through S 0 and L choices rather than from the crust densities themselves.

Systematic Error Estimate
In figure 3 (a) we plot the 2D Λ 1 vs Λ 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) 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.

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  figure 3 (b)). Right: Table containing the R X i radii values, obtained by matching different crusts to the X = 5 th or X = 95 th pressure-density percentile curves, inverting the implied M(R) relation and fixing M to M 1 = 1.57M or M 2 = 1.20M ( figure 4 (b)).  Figure 3: The Λ 1 vs Λ 2 (a) and mass vs radius (b) 2D distributions. The continuous and dashed curves represent, respectively, the 90% and 50% credible limits. While the Λ distributions are all indistinguishable -and the "prediction" curves perfectly superimposed! -the radii distributions obtained with the 30-65 and 34-35 crusts are systematically shifted with respect to each other and to the one resulting from the SLy low-density model. This shift measures the additional uncertainty, in radius only, due to the unknown crust EOS. The GW constraints on tidal deformation are insensitive to the crust.
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. This can be easily explained considering that for low-mass stars the crust has more relative importance than for heavier NS. Note also how our prediction (orange dashed) gives a good approximation of the curve obtained with the actual run (green)