Neutron Star Properties: Quantifying the Effect of the Crust–Core Matching Procedure

Abstract

The impact of the equation of state (EoS) crust-core matching procedure on neutron star (NS) properties is analyzed within a meta-modeling approach. Using a Taylor expansion to parametrize the core equation of state (EoS) and the SLy4 crust EoS, we create two distinct EoS datasets employing two matching procedures. Each EoS describes cold NS matter in a $\beta$ equilibrium that is thermodynamically stable and causal. It is shown that the crust-core matching procedure affects not only the crust-core transition but also the nuclear matter parameter space of the core EoS, and thus the most probable nuclear matter properties. An uncertainty of as much as 5% (8%) on the determination of low mass NS radii (tidal deformability) is attributed to the complete matching procedure, including the effect on core EoS. By restricting the analysis, imposing that the same set of core EoS is retained in both matching procedures, the uncertainty on the NS radius drops to 3.5% and below 1.5% for $1.9M_{\odot }$. Moreover, under these conditions, the crust-core matching procedure has a strong impact on the Love number $k_2$, of almost 20% for $1.0M_{\odot }$ stars and 7% for $1.9M_{\odot }$ stars, but it shows a very small impact on the tidal deformability $\mathrm{\Lambda }$, below 1%.

1. Introduction

Neutron stars (NSs) are astrophysical objects made of cold super-dense neutron-rich nuclear matter. They are unique systems through which one can access the properties of the equation of state (EoS) of nuclear matter under extreme conditions of density and isospin asymmetry. Some NS observations have imposed strong constraints on the EoS of nuclear matter, in particular, the masses of the pulsars PSR J1614−2230 with $M=1.908\pm 0.016 M_{\odot }$ [1,2,3], PSR J0348+0432 with $M=2.01\pm 0.04M_{\odot }$ [4], and the recently detected MSP J0740+6620 with $M=2.14{ }_{+0.10}^{-0.09} M_{\odot }$ [5]. The description of two solar mass pulsars requires a quite stiff EoS at large densities and may even inhibit the appearance of non-nucleonic degrees of freedom as discussed in e.g., [3,6]. However, many other works have shown that due to the rather limited number of constraints at high baryonic densities it is possible to predict NS with non-nucleonic degrees of freedom e.g., [7,8,9,10,11,12,13,14].

A strong constraint on the EoS would be set by the simutaneous measurement of the mass and radius of a NS. The already operating or programed missions NICER [15], Athena X-ray telescope [16], and the eXTP [17] aim precisely to measure the mass and radius of a NS with an uncertainty of the order of 5%. Recently, NICER has published their first results but the aimed accuracy was still not attained [18,19]. The gravitational waves (GWs) emitted during the coalescence of binary NS systems carry important information on the high density properties of the EoS. The GW170817 event has settled an upper bound of $\tilde{\mathrm{\Lambda }}\le 800$ on the effective tidal deformability of the binary $\tilde{\mathrm{\Lambda }}$ [20]. Under minimal assumptions about the nature of the compact objects, a follow up reanalysis of the GW170817 event gave $\tilde{\mathrm{\Lambda }}={300}_{-230}^{+420}$ [21]. Furthermore, the tidal deformability of a $1.4M_{\odot }$ NS was estimated to be $70<{\mathrm{\Lambda }}_{1.4M_{\odot }}<580$ [20]. However, this last prediction did not require that the EoS of the NS is able to account for two solar mass stars.

With the possibility that stricter observational constraints will be set in the near future, it is important to be able to set uncertainties on the quantities calculated from theoretical models. One possible origin of these uncertainties on the determination of the star radius stems from the way NS EoS is built. Neutron star matter should be described by an unified EoS [22,23,24], where all density regions are calculated from the same nuclear interaction, see also [25,26,27]. However, non-unified EoS built piecewisely using distinct models for different density regions [28,29,30,31] are frequently used. To construct a non-unified EoS, it is important to determine how and where the different density segments are matched. In fact, it has been shown that the matching procedure of the crust EoS to the core EoS may give rise to uncertainties as large as 10% [26].

In order to determine the constraints set by the gravitational waves GW170817 detected from the merging of two NS and other observations, several studies have been performed that apply a huge set of meta-model EoS [29,30,32,33,34,35]. In these studies different approaches have been considered to glue the core EoS to the crust. The most frequently used method was been applied in [28] to parametrize a well known set of EoS calculated from different theoretical approaches by a finite number of politropic functions. For all the models the crust determined in [36], based on the Skyrme interaction SLy4, has been considered and matched to the core at the pressure-energy density point where both EoS cross.

In the present work, we want to quantify how the matching of the crust EoS with the core affects NS properties, such as the radius and the tidal deformability. We will generate two datasets of non-unified EoSs that are characterized by different core-crust matching procedures: the first one is based on the $P\left({\mu }_c\right)$ relation and the second one on the $P\left({ϵ}_c\right)$, where ${\mu }_c$ and ${ϵ}_c$ are the core-crust chemical potential and the energy density at the transition, respectively. The first method uses a Maxwell construction, i.e., the transition occurs at constant baryonic chemical potential ${\mu }_c$, while the transition happens at a constant energy density, ${ϵ}_c$, in the second method. This other method was the one applied in [28] to build politropic parameterizations of well known EoSs. An analysis of the matching effect on the consistency of the causal and thermodynamic properties of the EoS was performed in [26], and it was shown that the matching on the pressure-energy density plane is thermodynamically inconsistent. However, it may still be justified if the unprecision introduced is within the uncertainty of the approach. Other methods have been considered such as an interpolation between some fixed baryonic density inside the inner crust and the saturation density [29], but we will only investigate the first two presented above.

In order to undertake the present study, we use a Taylor expansion around the saturation density of symmetric nuclear matter for the core EoS, while for the crust we choose the SLy4 EoS [36], since this was the crust EoS used in most of the works that applied meta-models.

The paper is organized as follows. In Section 2, we introduce the EoS parametrization and the crust-core matching procedures used in this work. The impact of crust-core matching procedure on several NS properties is analyzed in Section 3. The conclusions are drawn in Section 4.

2. EoS Parametrization

We start from the generic functional form for the energy per particle of homogeneous nuclear matter

$$\mathcal{E}(x,\delta )=e_{\mathrm{sat}}\left(x\right)+e_{\mathrm{sym}}\left(x\right){\delta }^2$$

(1)

with

$$e_{\mathrm{sat}}\left(x\right)=E_{\mathrm{sat}}+\frac{1}{2}{K}_{\mathrm{sat}}{x}^2+\frac{1}{6}{Q}_{\mathrm{sat}}{x}^3$$

(2)

$$e_{\mathrm{sym}}\left(x\right)=E_{\mathrm{sym}}+L_{\mathrm{sym}}x+\frac{1}{2}{K}_{\mathrm{sym}}{x}^2$$

(3)

where $x=(n-n_0)/\left(3n_0\right)$. The baryon density is given by $n=n_n+n_p$ and $\delta =(n_n-n_p)/n$ is the asymmetry, with $n_n$ and $n_p$ being the neutron and proton densities, respectively. This approach of Taylor expanding the energy functional up to fourth order around the saturation density, $n_0$, has been applied recently in several works [29,37,38,39].

The empirical parameters can be identified as the coefficients of the expansion. The isoscalar empirical parameters are defined as successive density derivatives of $e_{\mathrm{sat}}$,

$$P_{\mathrm{is}}^k={\left(3n_0\right)}^k{\left.\frac{{\partial }^k{e}_{\mathrm{sat}}}{\partial n^k}\right|}_{\{\delta =0,n=n_0\}},$$

(4)

whereas the isovector parameters measure density derivatives of $e_{\mathrm{sym}}$,

$$P_{\mathrm{iv}}^k={\left(3n_0\right)}^k{\left.\frac{{\partial }^k{e}_{\mathrm{sym}}}{\partial n^k}\right|}_{\{\delta =0,n=n_0\}}.$$

(5)

The corresponding empirical parameters are then

$$\{P_{\mathrm{is}}^0=E_{\mathrm{sat}},P_{\mathrm{is}}^2=K_{\mathrm{sat}},P_{\mathrm{is}}^3=Q_{\mathrm{sat}}\}$$

(6)

and

$$\{P_{\mathrm{iv}}^0=E_{\mathrm{sym}},P_{\mathrm{iv}}^1=L_{\mathrm{sym}},P_{\mathrm{iv}}^2=K_{\mathrm{sym}}\}$$

(7)

The saturation energy $E_{\mathrm{sat}}$ and saturation density $n_0$ being rather well constrained, we fix their values throughout this work: $E_{\mathrm{sat}}=-15.8$ MeV and $n_0=0.155$ fm${}^{-3}$.

Each possible EoS is represented by a point in the 8-dimensional space of parameters. We use random sampling of models through a multivariate Gaussian with zero covariance:

$$\begin{array}{cc}\hfill {\mathrm{EoS}}_i& ={(E_{\mathrm{sym}},L_{\mathrm{sym}},K_{\mathrm{sym}},K_{\mathrm{sat}},Q_{\mathrm{sat}})}_i\hfill \\ \hfill & \sim N(\mathit{\mu },\mathbf{\Sigma })\hfill \end{array}$$

(8)

where

$$\begin{array}{c}{\mathit{\mu }}^T=({\overline{E}}_{\mathrm{sym}},{\overline{L}}_{\mathrm{sym}},{\overline{K}}_{\mathrm{sym}},{\overline{K}}_{\mathrm{sat}},{\overline{Q}}_{\mathrm{sat}})\\ \mathbf{\Sigma }=diag({\sigma }_{E_{\mathrm{sym}}},{\sigma }_{L_{\mathrm{sym}}},{\sigma }_{K_{\mathrm{sym}}},{\sigma }_{K_{\mathrm{sat}}},{\sigma }_{Q_{\mathrm{sat}}}).\end{array}$$

(9)

In the present approach, as discussed in [29], no a priori correlations exist between the different parameters of the EoS. However, the requirement that every valid EoS must satisfy a set of experimental and observational constraints induces correlations among the parameters in the final set of EoS. The used parameters values and their standard deviations are listed in

Table 1. The mean ${\overline{P}}_i$ and standard deviation $\sqrt{{\sigma }_{P_i}}$ of the multivariate Gaussian, where ${\sigma }_{P_i}$ is the variance of the parameter $P_i$. Our equation of states (EoSs) are sampled using the initial distribution for $P_i$ assuming that there are no correlations among the parameters. All the quantities are in units of MeV.

${\mathit{P}}_{\mathit{i}}$	${\mathit{E}}_{\mathrm{sym}}$	${\mathit{L}}_{\mathrm{sym}}$	${\mathit{K}}_{\mathrm{sat}}$	${\mathit{K}}_{\mathrm{sym}}$	${\mathit{Q}}_{\mathrm{sat}}$
${\overline{P}}_i$	32	60	230	−100	300
$\sqrt{{\sigma }_{P_i}}$	2	15	20	100	400

.

We impose the following conditions to get a valid EoS: (i) the pressure is an increasing function of density (thermodynamic stability); (ii) the speed of sound is smaller than the speed of light (causality); (iii) the EoS supports a maximum mass at least as high as $1.97M_{\odot }$ [1,2,3,4] (observational constraint); and (iv) the symmetry energy $e_{\mathrm{sym}}\left(n\right)$ is positive. This may be a too restrictive constraint and a more realistic one would be that the symmetry energy $e_{\mathrm{sym}}\left(n\right)$ is positive for densities below the central density of the maximum mass star configuration. We consider, however, that the difference between both sets of EoS will not be significant. All EoS describe $npe\mu$ matter in $\beta$-equilibrium at zero temperature.

Crust Matching Procedure

We are going to compare two crust-core matching procedures used in the literature for obtaining the neutron star matter EoS. The first one consists in doing the matching in the $P\left(\mu \right)$ diagram, i.e., by requiring $P_{\mathrm{crust}}\left({\mu }_t\right)=P_{\mathrm{core}}\left({\mu }_t\right)$, where ${\mu }_t$ is transition baryonic chemical potential. The second method matches the crust and core EoSs in the $P\left(ϵ\right)$ diagram, i.e., by requiring $P_{\mathrm{crust}}\left({ϵ}_t\right)=P_{\mathrm{core}}\left({ϵ}_t\right)$, where ${ϵ}_t$ denotes de transition energy density. We further require that the crust-core transition occurs below $n_t<0.10$ fm${}^{-3}$, which is consistent with the range of core-crust transition densities for a large set of nuclear models [40].

3. Results

Using the same initial probability distribution function, we start by sampling an EoS as ${\mathrm{EoS}}_i\sim N(\mathit{\mu },\mathbf{\Sigma })$ and analyze whether: (i) the above conditions are met and (ii) a matching with the crust is possible. If both conditions are met, the ${\mathrm{EoS}}_i$ is taken as a valid EoS for $\beta –$equilibrated neutron star matter and added to the dataset.

3.1. Empirical Parameters Distribution

In this first section, the sampling procedure is independently done for each matching procedure, i.e., we are not testing if a sampled EoS gives a valid matching within both procedures. Instead, both matching procedures are tested on different and independently sampled EoSs. This way, we are able to study if the matching procedure affects the most probable set of parameters, i.e., the final distributions of the empirical parameters.

The final sets we have generated are made of 1326 and 1410 EoSs using the $P\left(\mu \right)$ and $P\left(ϵ\right)$ procedures, respectively. The parameters statistics of both sets are given in

Table 2. The mean and standard deviation (std) of the empirical parameters. All the quantities are in units of MeV.

	$\mathit{P}\left(\mathit{\mu }\right)$		$\mathit{P}\left(\mathit{ϵ}\right)$
	mean	std	mean	std
$K_{\mathrm{sat}}$	225.32	17.47	227.75	19.24
$Q_{\mathrm{sat}}$	−83.20	30.26	−83.69	34.93
$E_{\mathrm{sym}}$	33.18	1.84	31.41	1.90
$L_{\mathrm{sym}}$	61.89	8.14	77.70	7.47
$K_{\mathrm{sym}}$	−26.63	30.45	−32.32	37.07

. The mean values of the isoscalar properties are close, while the isovector show some differences. We show the density distribution of the isovectors properties in Figure 1. The main difference occurs for $L_{\mathrm{sym}}$, where the $P\left(\mu \right)$ procedure (blue) predicts lower values. It is already clear that the most probable EoS depends on the selected matching method. It is highly improbable that a specific core EoS will glue with a crust within both procedures.

The histograms for the transition densities, $n_t$, are shown in Figure 2. The values are concentrated at low densities for $P\left(\mu \right)$, while they spread over the whole density range for $P\left(ϵ\right)$ (additionally, we have imposed $n_t>0.04$ fm${}^{-3}$ as lower bound). In fact, if we soften the constraint $n_t<0.1$ fm${}^{-3}$ to $n_t<n_0=0.155$ fm${}^{-3}$, it would still result in similar low $n_t$ values for $P\left(\mu \right)$, while it would be highly clustered around $n_0$ for $P\left(ϵ\right)$. Therefore, both methods give consistently different transition densities for the crust-core transition. Moreover, under the $P\left(ϵ\right)$ matching procedure, the dataset statistics for the empirical parameters strongly depends on the constraint imposed on $n_t$.

In Figure 3, we show how the empirical parameters are correlated for each matching procedure. In the left panel, the upper triangular part indicates the correlation coefficient (the correlations are calculated via the Pearson coefficient $\mathrm{Corr} [x,y]=E\left[(x-{\mu }_x)(y-{\mu }_y)\right]/\left({\sigma }_x{\sigma }_y\right)$) for $P\left(ϵ\right)$, while the bottom part represents the $P\left(\mu \right)$. The correlation difference between both matching procedures is shown in the right panel. Is it clear that the matching procedure gives rise to substantial differences in the correlations between parameters, probably due to the differences occurring in the distribution of parameters. The major distinctions are seen in some of the stronger correlations, i.e., $\mathrm{Corr}[L_{\mathrm{sym}},K_{\mathrm{sat}}]$ and $\mathrm{Corr}[K_{\mathrm{sym}},K_{\mathrm{sat}}]$, while the smaller correlation difference happens for the weaker one, $\mathrm{Corr}[E_{\mathrm{sym}},K_{\mathrm{sat}}]$. Note, however, that the strongest correlation is $\mathrm{Corr}[L_{\mathrm{sym}},K_{\mathrm{sym}}]$ for both sets.

The predictions obtained from these two sets of models for the NS radius (left), Love number $k_2$ (center), and tidal deformability $\mathrm{\Lambda }$ are shown in Figure 4 as a function of the NS mass. The top panels show the mean (solid lines), two times the standard deviation (filled regions), and max/min values (dashed lines), where the color identifies the matching procedure: Blue for $P\left(\mu \right)$ and red for $P\left(ϵ\right)$. The middle panel shows the difference between means $\mathrm{\Delta }A=|{\overline{A}}_{ϵ}-{\overline{A}}_{\mu }|$, with ${\overline{A}}_{ϵ}$ and ${\overline{A}}_{\mu }$ denoting the mean values for the $P\left(ϵ\right)$ and $P\left(\mu \right)$ matching procedures sets, respectively. The bottom panels show the relative difference (%), defined by $\mathrm{\Delta }A=100\times |{\overline{A}}_{ϵ}-{\overline{A}}_{\mu }|/\overline{A}$, where $\overline{A}=({\overline{A}}_{ϵ}+{\overline{A}}_{\mu })/2$. We are analyzing how far is the average prediction of each set from each other.

The difference on the average predictions decreases as more massive NSs are considered. The average deviation between both predictions on the NS radius is of 5% for a $1.0M_{\odot }$ NS, and it decreases to <2 % for a $1.9M_{\odot }$ NS. For a $1.4M_{\odot }$ NS, the difference is smaller than $0.4$ km on average. The Love number predictions show the highest discrepancy: 18% and 6% for a $1.0M_{\odot }$ and $1.9M_{\odot }$ NS, respectively. The deviations on the $\mathrm{\Lambda }$ predictions can be as high as 8%. These deviations reflect the different $P\left(ϵ\right)$ regions span by each set, or, more precisely, the gap between the most probable regions, since the matching procedure affects the core EoS. Therefore, when a similar statistical approach is used for generating a core EoS, the matching procedure choice carries the above uncertainties on the NS properties.

The effective tidal deformability of a binary system is the leading tidal parameter of the gravitational-wave signal from a NS merger. It is given by

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

(10)

where $q=M_2/M_1<1$ is the binary mass ratio and ${\mathrm{\Lambda }}_1 \left(M_1\right)$ and ${\mathrm{\Lambda }}_2 \left(M2\right)$ represent the tidal deformability (mass) of the primary and the secondary NS in the binary, respectively. The GW170817 event provides an upper bound of $\tilde{\mathrm{\Lambda }}={300}_{-230}^{+420}$ (low spin-prior) [21], while the binary mass ratio is bounded as $0.73\le q\le 1$ [21]. The binary chirp mass,

$$M_{\mathrm{chirp}}=\frac{{\left(M_1{M}_2\right)}^{3/5}}{{(M_1+M_2)}^{1/5}}=M_1\frac{q^{3/5}}{{(1+q)}^{1/5}},$$

(11)

is a another quantity that is measured with a great accuracy from the gravitational-wave signal. It was measured to be $M_{\mathrm{chirp}}=1.{186}_{-0.001}^{+0.001}{M}_{\odot }$ [21] for the GW170817 event.

In Figure 5, we have determined the effective tidal deformability, $\tilde{\mathrm{\Lambda }}$, of all binary systems compatible with $M_{\mathrm{chirp}}=1.186M_{\odot }$ and $0.7<q<1$. Both sets of EoSs show very similar results: A small weak dependence of $\tilde{\mathrm{\Lambda }}$ on q (as already noted in [30,41]), and the deviations due to the crust matching procedure, as measured by the difference on their mean value predictions, are below 6%. This result indicates that when a dataset of EoS is constructed, from a nuclear EoS modeling or a meta-modeling approach, with a specific crust-core matching procedure, an uncertainty as high as 6% must be taken into account in the $\tilde{\mathrm{\Lambda }}$ predictions.

Let us stress that the uncertainties discussed above result from two different contributions: (i) the crust-core matching procedure itself and (ii) the implication that the crust-core matching procedure has on the allowed parametrization space of the core EoS (see Table 2 for the differences). The two contributions are entangled in the present analysis. In the next section we will consider the first one alone.

3.2. Isolating the Matching Procedure Effect

In this section, we isolate the impact of the matching procedure on the NS properties. For that, we are interested in finding a set of EoS that simultaneously allows for a crust-core transition within both matching procedures. From the previous section, we have learned, however, that it is quite unlikely that a specific EoS sampled from the initial probability distribution would enable a crust-core matching within both procedures, as the considerable difference on the parameters distribution and correlations seems to indicate. In other words, the overlap region between the two final probability distributions is very small, indicating that, once an hadronic EoS is chosen, it is unlikely that both matching procedures are possible. However, by drawing a considerable number of EoS samples from the initial probability distribution, we were able to construct a set of 55 EoS that satisfy both matching procedures.

In Figure 6, we show the mean, and some other statistics, of the prediction difference for NS radius (left), Love number $k_2$ (center), and tidal deformability $\mathrm{\Lambda }$. Let us point out that in Figure 4 we have defined $\mathrm{\Delta }R$ (e.g.) as the absolute difference between means, $|{\overline{R}}_{ϵ}-{\overline{R}}_{\mu }|$, and in Figure 6, we are considering $\mathrm{\Delta }R=\frac{1}{N}{\sum }_i|R_{ϵ}^{{\mathrm{EoS}}_i}-R_{\mu }^{{\mathrm{EoS}}_i}|$ for the prediction difference and $100\times (\frac{1}{N}{\sum }_i|R_{ϵ}^{{\mathrm{EoS}}_i}-R_{\mu }^{{\mathrm{EoS}}_i}\left|\right)/\overline{R}$, where $\overline{R}=({\overline{R}}_{ϵ}+{\overline{R}}_{\mu })/2$, for the relative difference. Thus, $\mathrm{\Delta }R$ denotes here the average radii difference instead of the difference on the radii averages as in the last section.

We see that the discrepancy on the NS radius prediction lies between 2.5% and 5% corresponding to ∼0.34 km and ∼0.68 km for $1.0M_{\odot }$, and it decreases to 1% and 2% corresponding to 0.1 km and 0.25 km for $1.9M_{\odot }$. The Love number $k_2$ is the quantity that shows the highest prediction difference, with almost 20 for $1.0M_{\odot }$ and 7% for $1.9M_{\odot }$. On the other hand, regardless of the NS mass, the crust-core matching method has only a small impact on the $\mathrm{\Lambda }$, lower than 1%. Considering that $\mathrm{\Lambda }=\frac{2}{3}{k}_2{C}^{-5}$, where $C=M/R$ is the star compactness, the quite small $\mathrm{\Delta }\mathrm{\Lambda }$ value is surprising when compared with both $\mathrm{\Delta }R$ and $\mathrm{\Delta }{k}_2$. This result rises from a cancellation effect between $k_2$ and C due to the non-trivial dependence of $k_2$ on C and on the EoS [27,42,43]. Despite the sensitivity of $k_2$ and R on the matching procedure for the same core EoS, the tidal deformability shows no such degree of sensitivity. From this analysis, when considering a non-unified EoS approach in describing neutron star matter, we are able to attribute a prediction error that arises solely from crust-core matching procedure.

Finally, in Figure 7, we display the prediction difference for $\tilde{\mathrm{\Lambda }}$ as a function of q (fixing $M_{\mathrm{chirp}}=1.186M_{\odot }$ and $0.7<q<1$). The crust-core matching procedure has only a marginal impact on the effective tidal deformability of NS binaries, which are compatible with the GW170817 event.

We conclude that NS properties predicted from a non-unified EoS approach, where a crust is added under a certain matching procedure, carry a systematic uncertainty, which includes the influence of the matching procedure on the accepted core EoS. An uncertainty of ∼5% was determined in the present study for the radius of NSs with masses below 1.4 $M_{\odot }$.

4. Conclusions

We have analyzed the impact of the EoS crust-core matching procedure on the NS properties. We have generated two distinct EoS datasets employing two matching procedures used in the literature, where the matchings are done in the $P\left(\mu \right)$ and $P\left(ϵ\right)$ planes. The core EoS was determined from a Taylor expansion approach, in which the energy functional is expanded in a Taylor series around the saturation density. Different core EoSs were generated through random sampling of the empirical parameters via a multivariate Gaussian with zero covariance. To get a valid EoS, we then impose thermodynamic stability, causality, and a maximum star mass of 1.97 $M_{\odot }$ to each sampled EoS.

In the first part, we have concluded that, when meta-modeling is used to generate a set of all possible nuclear matter EoS, the crust-core matching procedure affects the output EoS parameter space, and thus the NS properties will show distinct statistics concerning their properties. When considering the difference between the mean values of each dataset, we saw that the impact of the crust-core matching procedure on the NS properties decreases as more massive NSs are considered: an uncertainty of ∼5% was obtained for a $1.0M_{\odot }$ and below 2% for a $1.9M_{\odot }$ NS. The Love number predictions show a discrepancy of 18% and 6% for $1.0M_{\odot }$ and $1.9M_{\odot }$ NS, respectively. The deviations on the $\mathrm{\Lambda }$ predictions can be as high as 8% for a $1.0M_{\odot }$ NS but go down to a value below 4% for $M\ge 1.5M_{\odot }$. We notice that the uncertainty introduced in the calculation of the radius for stars with a mass $M\le 1.4M_{\odot }$ is of the order of the uncertainty that NICER aims to attain in the future observations.

In the second part, we isolate the impact of the matching procedure by building a dataset of core EoSs that allows for both matching procedures under the applied constraints, i.e., $0.04<n_t<0.1$ fm${}^{-3}$. We concluded that uncertainty in the NS radius prediction is in between 0.34 km and 0.68 km for $1.0M_{\odot }$, and it decreases for 0.1 km and 0.25 km for $1.9M_{\odot }$. The Love number $k_2$ shows the highest prediction difference, with almost 20% for $1.0M_{\odot }$ and 7% for $1.9M_{\odot }$. However, the crust-core matching method shows a very small impact on the $\mathrm{\Lambda }$ (below 1%). This conclusion is consistent with the results obtained in [27,42,43], when the analysis of the effect of the crust is done taking a fixed core EoS. Likewise, the effective tidal deformability values of the of NS binaries, compatible with the GW170817 event, show almost no effect from the crust-core matching procedure.

Let us mention that the EoS parameters in Equations (2) and (3) are seen as unknown parameters to be extracted from astrophysical observations, and not the Taylor coefficients from any energy density functional of nuclear matter [44]. We interpret our EoS as metamodels, and while the nuclear matter parameters of the lower order expansion terms may have a direct relation with the nuclear matter properties constrained by experimental measurements, the parameters connected to the higher order terms should be considered as effective parameters that will be constrained by the astrophysical observations [37,44].