{'fraction_non_alphanumeric': 0.03463005579899751, 'fraction_numerical': 0.019419312127612623, 'mean_word_length': 4.524949925977532, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Polytropic fits of modern and unified equations of state

Equations of state for a cold neutron star's interior are presented in three-column tables that relate the baryonic density, the energy density, and the pressure. A few analytical expressions for those tables have been established these past two decades, as a convenient way to present a large number of nuclear models for neutron star matter. Some of those analytical representations are based on nonunified equations of state, in the sense that the high and the low density part of the star are not computed with the same nuclear model. Fits of equations of state based on a piecewise polytropic representation are revised by using unified tables of equations of state, that is to say models which have been calculated consistently for the core and the crust. A set of 52 unified equations of state is chosen. Each one is divided in seven polytropes via an adaptive segmentation, and two parameters per polytrope are fitted to the tabulated equation of state. The total mass, radius, tidal deformability and moment of inertia of neutron stars are modelled from the fits and compared with the quantities calculated from the original tables to ensure the accuracy of the fits on macroscopic parameters. We provide the polytropes parameters for 15 nucleonic relativistic mean field models, seven hyperonic relativistic mean field models, five hybrid relativistic mean-field models, 24 nucleonic Skyrme models, and one ab initio model. The fit error on the macroscopic parameters of neutron stars is small and well within the estimated measurement accuracy from current and next generation telescopes.


I. INTRODUCTION
The equation of state (EoS) plays a key role in modeling a neutron star's macroscopic parameters. Despite continuous efforts to push the limits of nuclear experiments, conditions of density and temperature in the deepest layers of neutron stars remain out of reach for laboratories. Multimessenger astronomy provides a chance to probe deep inside these extremely compact stars: the relativistic hydrodynamics equations operate as a bridge between the unknown microphysics of highdensity neutron-rich matter, and observable macroscopic parameters. The measurement of the mass M , the radius R, the tidal deformability Λ, and the moment of inertia I paired with hydrodynamics equations in gravity theories, have turned neutron stars (NSs) into extra-terrestrial laboratories for high density nuclear physics.
Neutron stars are observed in various wavelengths of the luminosity spectrum. On one hand, they are observed in the x-ray spectrum. The spatial telescope X-ray Multi Mirror-Newton (XMM-Newton [1]) from the European Spatial Agency, and its American counterpart Chandra [2], have been operating for more than 20 years, providing data for isolated, accreting, and highly magnetized neutron stars. The younger Nuclear Spectroscopic Telescope Array (NuSTAR [3]), and Neutron-star Interior Composition ExploRer (NICER [4]) observe respectively in hard and soft x-rays. The next generation of x-ray telescopes is already on its way: the Enhanced X-ray Timing and Polarimetry (eXTP [5]) mission is scheduled to be launched in 2027, and the highly anticipated Advanced Telescope for High ENergy Astrophysics (ATHENA [6]) with a state-of-the-art X-ray Integral Field Unit (X-IFU) spectrometer will provide unprecedented spectral resolution for a wider effective area, hopefully in the 2030s. On the other hand, neutron stars are observed as radio sources, for example, in the observatories of Parkes, Green Bank, and Nancay. After almost 60 years of operation in Porto Rico, Arecibo fell in November 2020; fortunately, the largest worldwide telescope, the Square Kilometer Array (SKA [7]) will be operational in a few years and will include the 2017 launched Chinese contribution Five-hundred-meter Aperture Spherical radio Telescope (FAST [8]). Gravitational wave (GW) detection is the most recently explored area for compact object messaging: the LIGO Scientific Collaboration, the Virgo Collaboration, and the KAGRA Collaboration (LVK [9]) have provided promising results with the detection of doubleneutron-star binary mergers, and will keep on with run O4 starting in early 2023, and then O5. The Einstein 9 10 11 12 13 14 15 R (km) 1 Mass-radius (MR) sequence for the unified nucleonic and hyperonic relativistic mean-field EoSs used in this paper. Mass measurements of J0740+6620 and J1614−2230 within 1σ precision are presented, as well as NICER mass-radius contours of J0740+6620 and J0030+0451 within 1σ (see text for details), and the mass-radius constraint on GW170817. The exclusion regions for finite pressure and subluminal EoSs are represented respectively in yellow and green.
telescope [10], and its American counterpart the Cosmic Explorer [11], are both scheduled to start observing in the mid 2030s, and should detect continuous gravitational waves.
Part of the outer crust of neutron stars is well constrained by laboratory measurements of mildly neutron rich nuclei. However, the inner crust and the core are respectively poorly, and very poorly constrained which leads to many different supra-nuclear theories: nucleonic matter, confined or deconfined quarks, hyperonic matter, strange quark matter, etc. There exists dozens of equations of state for cold and highly dense neutron-star matter, a large number of them are gathered in an online data base named CompOSE [12]. The crust equation of state is more difficult to compute than the core's because it requires a treatment of inhomogeneous matter. For that reason, one can find nonunified constructions of equations of state: the core and the crust are not computed using the same nuclear model; the common practice is to attach an already established crust to the core equation of state. Nonunified constructions are oftentimes found in the literature, but have been shown to result in large errors in the modeling of macroscopic parameters, see Ref. [13].
The output format for computations of cold matter equations of state is a three column table with the baryonic density n, the energy density ǫ, and the pressure P .
However, an analytical form of tabulated equations of state is convenient, particularly for neutron-star simulations. To establish analytical representations of equations of state, one chooses a parametrized expression, and then adjusts parameters to the tabulated equation of state. Having one expression with easily comparable parameters is also a practical way to compare microscopic and macroscopic features of neutron stars. Adjusting the parametrized representation is called a fit: a few have been proposed so far, e.g., spectral fits as presented in Ref. [14], and piecewise polytropic fits by Ref. [15].
In this paper, piecewise polytropic fits based on nonunified constructions are revised, and new fits for 52 unified equations of state of neutron star's matter are presented. In Sec. II, the most relevant microscopic and macroscopic features required for a realistic equation of state are presented. After a brief overview of the various frameworks for nuclear interaction modeling used to compute equations of state, we present the equations of state used in this paper. Piecewise polytropic fits as well as the relation between polytropic parameters, and quantities n, ǫ and P , are presented in Sec. III. The importance of using unified equations of state, on the accuracy of macroscopic parameters is investigated; this section ends with details on our fitting method. In Section. IV, the accuracy of fits on the mass, the radius, the moment of inertia, and the dimensionless tidal deformability calculated within the framework of Einstein's theory of general relativity is presented for key astrophysical quantities. The role of nonunified constructions on the accuracy of socalled "universal" relations is discussed. The parameters for the fits are presented in tables in Appendix V.

II. CHOICE OF EQUATIONS OF STATE
The three basic physics rules that an EoS is required to meet are: • thermodynamic consistency: the first law of thermodynamics must be fulfilled, • causality: interpreted as subluminal EoS (v sound < c), for details concerning Lorentz invariance and causality see Ref. [16], • Le Chatelier's principle which states that the energy increases with pressure.
The natural limit up to which one can study the properties of dense matter using NS observations is set by the properties of matter at the center of the star with maximum mass configuration -the corresponding value of the central baryon number density is denoted n max . Astrophysical and nuclear laboratory data can impose more constraints to ensure that the equation of state is realistic for NS matter. The macroscopic parameters of cold NSs are computed by injecting the EoS into the equations of hydrodynamics equilibrium within a given theory of gravity. Here we adopt Einstein's theory of general relativity. Comparing modelled macroscopic parameters to their measurement provides information on the NS matter EoS. More to the point, the quality of an EoS relies on its ability to be consistent with observations of the NS's astrophysical quantities.
The modeling of a NS's macroscopic parameters treated within general relativity imposes a first constraint on the mass/radius ratio: it must always be larger than that of Schwarzschild which excludes any mass/radius relation for which R < 2GM/c 2 . A second constraint is imposed by the finite nature of the pressure, and renders R > 9GM/4c 2 (yellow region in Figs. 1,2). A third constraint on the mass/radius ratio is imposed to ensure a subluminal EoS: R > 2.9GM/c 2 (green region in Figs. 1,2,3).
An important test for an EoS is to reach the highest NS mass measured: in this article, we use the massive millisecond pulsar J1614−2230 [17] measured with a mass of 1.908 ± 0.016 M ⊙ (solar mass M ⊙ ) as a selection criterion. Sources with larger masses have been reported, such as the millisecond pulsar J0740+6620, whose mass was previously measured with relativistic Shapiro delay at 2.14 +0. 10 −0.09 M ⊙ [18], but was recently revised to a 2.08±0.07 M ⊙ [19]. The source J0348+0432 was reported with a mass of 2.01 ± 0.04 M ⊙ [20] and J1810+1744 has been measured with a mass of 2.13 ± 0.04 M ⊙ [21]; their mass measurement technique is based on a highly model dependent analysis of the companion's photometry (white dwarf). We refer to Fig. 1 of Ref. [13] for an overview of NS measured masses. The maximum mass of a cold nonrotating star was constrained by gravitational wave (GW) detection [22] of the double NS binary GW170817 in Ref. [23], which highlighted the importance of finite temperature in the relation between the gravitational mass Keplerian limit and the maximum mass of a nonrotating cold star.
The radius is directly connected to the EoS's stiffness. For a significant range of pressures typical of neutron-star interior, the corresponding densities are smaller for a stiff EoS than for a soft one. As a result, from the Tolman-Oppenheimer-Volkoff equations, the thickness of this region is larger for stiff EoSs which can lead to a larger radius of the NS. A strong correlation exists between the pressure and the radius at densities [1 − 2.5]n 0 (saturation density n 0 = 0.16 fm −3 ), as shown by Ref. [24]. The relations between the total mass and the total radius for the EoSs used in this paper are presented in Figs. 1, 2 and 3. The radius and mass of two sources -PSR J0030+0451 and PSR J0740+6620, have been reported by the NICER telescope by two teams each. The measurement technique is based on an analysis of the surface emission of the pulsar, precisely of its hot spots. The source J0030+0451 was reported by Ref. [25] to have a mass of 1.34 +0. 15 −0.16 M ⊙ and radius of 12.71 +1.14 −1.19 km and reported by Ref. [26] to have a mass of 1.44 +0. 15 −0.14 M ⊙ and a radius of 13.02 +1.24 −1.06 km within 1σ precision. The source J0740+6620 was reported by [27] to have a radius of 13.7 +2.6 −1.5 km and by [28] to have a radius of 12.39 +1.3 −0.98 km within 1σ precision; prior knowledge from XMM-Newton telescope on the mass for this source was used. Contours for those sources are presented in Figs. 1, 2 and 3; for the source reported by Ref. [26], the contours for two signal analysis are given because no preference for one or the other was significant. Because the uncertainty for the radius is quite large, measurements serve more as a proof of concept for an elegant radius determination, than a conclusive constraint on NS matter. In Refs. [29,30], an indirect estimation of the radius was established from GW170817 gravitational wave (GW) detection; the authors either used "universal" relations established by Ref. [31] to present a radius with 3.5 km error bars, or a collection of EoSs in a Bayesian analysis to give a likelihood for the radius. As a substitute for radius measurement, a series of papers attempt to impose limits on the radius. Reference [32] used a prior distribution of EoSs to obtain a radius interval of [10.4−12.9] km for a 1.4 M ⊙ NS with a 95% confidence level. A similar approach is used by Ref. [33] in a Monte Carlo analysis with five low mass x-ray binary sources, to extract a minimum radius of 9.1 +1.3 −1.5 km within a 90% confidence level. Unfortunately, those limits are established from a small set of sources. Reference [34] established a constraint on the radius, based on the assumption that NS's rotation follows a Keplerian frequency; the authors conclude that R 1.4 ≤ 9.3 km for the radius of a 1.4 M ⊙ . However, NSs might not follow such frequencies when the rotation of the star is associated with GW emission or when triaxial deformability sets on.
The tidal deformability Λ is the propensity of a star to be deformed by a neighboring gravitational field. For the binary NS merger GW170817, the tidal deformabilities of the stars were extracted from the inspiral waveform. This detection has indicated a preference towards a soft EoS, see [35,36]. Conciliating a high enough maximum mass with a soft EoS, while keeping a small radius is a delicate question of balance for the inner core composition. Constraints on the tidal deformability of a 1.4 M ⊙ NS are discussed in Ref. [29].
The moment of inertia I is a parameter that has not been measured so far. It requires a monitoring of the relativistic features of the binary orbit (the more compact the better) over a long period of time, see [37]. This parameter would be best measured in a pulsar binary such as the famous PSR J0737−3039 [38]. One can extract I using "universal" relations: [39] have used a Markov chain Monte Carlo analysis of NICER J0030+0451 data to establish a radius distribution for the star of mass 1.3381 M ⊙ in double binary PSR J0737−3039; from the compactness C = GM/Rc 2 , they estimate the moment of inertia.

B. What makes a reasonable NS EoS from the microscopic point of view ?
A large portion of the outer crust is constrained by laboratory measurements of mildly rich nuclei; for example, Ref. [40] has catalogued such results in a table with thousands of nuclei data. All models for NS interior should therefore be similar for the first few layers of the crust.
The symmetry energy and the slope of the symmetry energy at saturation density, respectively J and L, are two microscopic parameters used to describe an EoS. The symmetry energy is defined as the difference between the energy per baryon calculated for pure neutron matter, and the energy per baryon calculated for symmetric matter. A series of laboratory experiments can constrain both those parameters; we refer to Table II and Fig. 8 of Ref. [41] for a systematic approach. We refer to Ref. [42] for a compilation of constraints on the symmetry energy obtained from experiments and theory. In the following, we briefly present, nonexhaustively, how both these parameters can be constrained by laboratory data.
The binding energy described by the Finite Range Droplet Model (FRDM) includes symmetry-related terms whose values can be explored by using large ta- MR sequence for the unified hybrid relativistic mean-field EoSs used in this paper. Mass measurements of J0740+6620 and J1614−2230 within 1σ precision are presented, as well as NICER mass-radius contours of J0740+6620 and J0030+0451 within 1σ (for details see text), and the mass-radius constraint on GW170817. NL3 models are in blue and DD2 models are in green. The pressure bag constant B is either set to zero (plain lines) or takes a nonzero value (dotted lines). The exclusion region for subluminal EoSs is represented in green.
bles of nuclei data. The formula for the FRDM binding energy presented in Ref. [43] is paired with the table of nuclei of Ref. [44] to extract: L = 70 ± 15MeV and J = 32.5 ± 0.5MeV. Coulomb effects are linked to the surface symmetry term in the FRDM (see Ref. [45], Sec. IIC); experiments are performed on isobaric nuclei [46] in order to alleviate this entanglement. The same method has been used with Skyrme forces in Ref. [47]. The Heavy Ion Collision (HIC) of nuclei such as gold has introduced constraints on symmetric matter over the saturation density. The collision of nuclei such as isotopes of tin allows one to probe the asymmetry between the number of protons and the number of neutrons [48]. Recently, the spectral pion ratio of tin isotopes has been used to determine the slope L of the symmetry energy at saturation density and, at a 95% confidence level, 42 < L < 117 MeV, an interval consistent with the conclusions drawn in Ref. [42].
Neutron-rich nuclei are particularly interesting to investigate the symmetry energy when they present an asymmetric number of neutrons and protons, as is the case of Sn isotopes or 208 Pb that closes its nucleon shells, which simplifies the nucleus structure. An asymmetry in favor of neutrons implies that the nucleus will present a large difference in the radius distribution of neutrons and

Oertel+17
FIG. 4: Symmetry energy and its slope at saturation density respectively denoted J and L for Skyrme models used in this paper. Experimental data constraints are presented: in blue is the compiled constraint presented in Ref. [41], and in red that of PREX-II. The names of the EoSs in green refer to nucleonic models that do not permit the Direct Urca process.
protons, also called neutron skin. There are a few different ways of measuring the neutron skin thickness, one of which is to see how electroweak parity of 208 Pb is violated by polarized electrons in the experiments PREX-I and PREX-II [49,50] and of 48 Ca in the experiment CREX [51]. To extract L from neutron skin thickness measurements, the correlation between the two quantities is exploited via a fit established within a nuclear theoretical framework. In Ref. [52], the Skyrme Hartree-Fock model is used on measurements of Sn isotopes to constrain the relation between J and L. In Ref. [53], the FSU2Gold relativistic mean-field parametrization is used on PREX-II data to extract J = 38.1 ± 4.7 MeV and L = 106 ± 37 MeV; this result is in tension with other nuclear experiment constraints. An analysis of the compatibility between PREX-I and PREX-II and CREX experiments and other experiments determining J and L is discussed in Ref. [54]. The values of J and L for Skyrme and relativistic mean-field models used in this paper are presented respectively in Figs. 4 and 5. The laboratory constraints from Ref. [41], and PREX-II data are shown. Constraints on L and J can be combined with NS luminosity observations to explore the core composition. Indeed, Direct Urca (DUrca) which is a neutrino-emissive rapidly cooling process, is required to explain the cooling of accreting NSs [55]. This process is permitted if the proton fraction is high enough, and is therefore

PREX-I I Oertel+17
FIG. 5: Symmetry energy and its slope at saturation density respectively denoted J and L for relativistic mean-field models used in this paper. Experimental data constraints are presented: in blue is the compiled constraint presented in [41], and in red that of PREX-II. EoS's name in green refer to nucleonic models that do not permit the Direct Urca process.
triggered at a given value of the density n DUrca . This threshold is constrained by the density dependence of the symmetry energy, such that a large L favors a large proton fraction and a process allowed for lower n DUrca -equivalently NS mass. The presence of hyperons in the core implies that there is no need for an elevated L to trigger the DUrca process, it necessarily appears for a n DUrca < n max , for details see e.g., Ref. [56]. In the set of EoSs presented in this paper, nucleonic models DD2, DDME2, BSk20, BSk26, KDE0v1, SLy2, DH and Skb do not permit DUrca (n DUrca > n max ); they are presented in green in Fig. 5 and 4. The collective motion of nuclei is a source of giant resonances: let there be an exterior isoscalar monopole operator, the strength function of excited states in response to that operator is directly linked to the nuclear incompressibility K for which experimental data are available, (see Tables I and II of Ref. [57]); again, the relation between experimental data and K is established within a theoretical framework (e.g., Skyrme or Gogny forces). Constraints on L and J can also be extracted, such as presented in Refs. [58,59].
There are also some attempts at including results from cold atom experiments to constrain a part of the low density EoS using the unitary Fermi gas approach. The idea is to consider that low-density neutron matter can be solely characterized by an infinite range scattering and can be considered as a unitary Fermi gas -see Chap. 2 of Ref. [60] for details-in which the energy per nucleon is determined by a single and universal parameter. There are, however, no considerations of lattice nor clusters with this approach, which are essential in the understanding of crust physics.
Techniques to infer knowledge on microscopic parameters of EoS from astrophysical measurements were recently designed. Bayesian inference statistical analysis of GW signal data helps provide constraints on the EoS; we refer to Ref. [61] for a recent study of this technique, and to Ref. [62] and reference therein for a review. GW170817 detection was used in Ref. [63] to constrain L: the effective tidal deformability of the binary is strongly correlated to that quantity and GW170817 indicates a preference for a low L. Reference [35] discovered strong correlations between the tidal deformability and linear combination of pairs of nuclear parameters of different orders. References [27,28] related to NICER observation of J0740+6620, used a Bayesian analysis based on likelihood of measurement to infer high density EoS properties. Perfect knowledge of the EoS under half the saturation density (n 0 /2 = 0.08 fm −3 ) is assumed; over n 0 /2, a parametrized EoS and Gaussian-process-based models are used. A similar technique is used by Ref. [64] with Chiral Effective Theory (CET) constraints. Note that these techniques depend strongly on the priors chosen.

C. The variety of core composition
In the innermost parts of NSs, densities can go as high as 15n 0 . Such a supranuclear framework is out of reach for laboratories which leaves the description of NS cores open to various composition hypotheses. In the present paper, we use the following three categories: nucleonic, hyperonic, and hybrid models for the core composition.
In addition to nucleons in the core, the presence of hyperons -baryons with at least one strange quark-softens the core EoS, which induces a smaller radius; hyperons in NSs were first introduced in the 1960s. Because the nucleonic Fermi pressure is higher than the hyperonic pressure at fixed density, allowing the presence of hyperons results in a smaller radius. Softening the EoS, however, leads to a lower maximum mass. This is a problem referred to as the "hyperon puzzle", which can be counteracted if one finds a way to stimulate hyperonic pressure. One way to do so is to instigate repulsion from the nature of baryon interactions, see Refs. [65,66]. Hyperons are expected to appear over 2n 0 through a series of reactions involving nucleons such as (but not exclusively) Those reactions have their inverse thus creating loops of high-neutrino-emissive processes; the role of hyperonic neutrino emission on the cooling of NSs was studied by Ref. [67]. Hyperons Λ(uds), Σ 0 (uds), and Σ − (dds) once created in Eqs. (1), can themselves be sources of double strange Ξ − (dss) hyperons. Laboratory measurements of hyperons are performed in heavy-ion collision experiments at the Thomas Jefferson National Accelerator facility (Jlab, USA), the Mainz Microtron accelerator (MAMI-C, Germany), and the Japan Proton Accelerator Research Complex (J-Parc, Japan); for details see Ref. [68]. However, only very short lived hypernuclei can be measured, making it difficult to give solid constraints on the parametrization of YN (hyperon-nucleon) and YY (hyperon-hyperon) interactions. For example, no scattering has been measured, which is needed to accurately calculate hyperon properties. Attempts at extracting the strange composition of NSs from observations of mass, radius and thermal evolution have concluded that the radius of the star decreases linearly with the increase of the total hyperon content, see Ref. [69].
In the inner core, a phase transition from confined (baryons) to deconfined quarks can be included: the core is therefore made of (strange-)quark matter and EoSs are called hybrid.
A common approach to compute hybrid EoSs is the relativistic mean-field (RMF) approach. Quark matter is frequently described in the framework of the Nambu-Jona-Lasinio (NJL) model, the quark-meson model or the MIT bag model. The hybrid EoSs we consider in the present study were built in Refs. [70,71] within the SU(3) NJL model for quark matter and a RMF model for hadronic matter. The confined phase follows the same Lagrangian density as nucleonic or hyperonic (strange quark) matter. The Lagrangian density of deconfined quarks in the NJL model includes four-quark scalar and pseudoscalar interaction terms with coupling constant G S , four-quark vector and pseudovector interaction terms -both vector-isoscalar (VP) with coupling constant G ω and vector-isovector (VIPI) with coupling constant G ρ are considered-, and the six-quark t'Hooft term that ensures that the axial symmetry U(1) A is broken. In the NJL models, the pressure and energy density are defined up to a constant B. The pressure bag constant B is chosen to ensure either that the effective pressure falls to zero when the baryon chemical potential vanishes (B0) or to impose another type of constraint such as fixing the deconfinement baryonic density (B). Finally, the ratio between the vector and scalar coupling constants, in particular, ξ = G ω /G s and η = G ρ /G s , are parameters that characterize the models and define the intensity of the VP and VIPI channels, respectively (for more details, see Refs. [70] and [72]). In the following, we designate the hybrid EOS by the hadronic EoS-Bx − 100ξ − 100η, where x characterizes the magnitude of the bag constant in MeV/fm 3 . The quark phase transition induces a density jump that ensures a softening of the EoS. If this jump exceeds a critical value, the softening is such that the MR sequence presents a branch which is partly unstable (with respect to radial oscillations). In this case, a single EoS corresponds to the two stable branches of stellar models, and twin stars with the same mass but different radius can exist.
Hybrid stars have also been described within ab-initio approaches such as the Brueckner-Hartree-Fock manybody theory with realistic two-body and three-body forces [73,74] (for other ab-initio approaches see the reviews [41,75]).

D. Set of equations of state
The large number of available EoSs stems from two variables: the core composition, and the nuclear theory used for computations. There are two approaches to compute particle interaction of NS matter. The microscopic one is based on many-body interactions: such models are time-consuming and only accurately apply to homogeneous matter. The phenomenological approach consists of theories based on parameters adjusted to reproduce the properties of matter measured in laboratories.
The RMF description of hadronic matter is a phenomenological approach in which particles are considered to be immersed in a common potential established in quantum field theory. The Lagrangian density for nucleonic matter contains terms for (naked) baryons, leptons and quarks, as well as terms for interacting mesons of type σ (scalar), ω (vector-isoscalar), and ρ (vectorisovector), with self and cross interactions, see Ref. [76].
The EoS is obtained from the Lagrangian density where the first term L i is the nucleonic term, the last term L lep is the leptonic contribution, and the middle terms refer to the mesonic contributions. The nucleon Lagrangian density reads with the covariant derivative iD µ = i∂ µ − g ω ω µ − gρ 2 τ · ρ µ , g σ , g ω and g ρ being the meson-nucleon couplings. τ the SU(2) isospin matrices, and the effective mass is M * = M − g σ σ , with M being the vacuum nucleon mass. The leptonic Lagrangian density is given by where the sum is over electrons and muons, and m i is their mass. The mesonic contributions are: , m σ , m ω , and m ρ are the meson masses and κ, λ, ζ, ξ and Λ ω are constant coupling parameters. We will consider models with constant couplings, NL, NL-ωρ, TM, TM-ωρ, GM, H, FSU; the BSR models which include nonlinear meson terms, and DD models which have density dependent couplings, and for which the couplings related to nonlinear mesons terms are zero (in particular, the couplings κ, λ, ξ, ζ and Λ ω ). For the DD models the isoscalar couplings of the mesons i to the baryons take the form and the isovector ones are given by In these expressions, n 0 is the symmetric nuclear saturation density, and x = n/n 0 . For more information, we refer to references listed in the last column of Table I. For hyperonic matter, terms for naked hyperons and mesons mediating the hyperonic interaction are added to the Lagrangian density. Coupling constants for the YY and YN interactions are presented in Refs. [69,76,77]. For the hyperonic sector, interactions are much simpler and include no cross terms between the hyperonic mesons [76]. For quark matter, terms for scalar, vector, and pseudovector quark couplings are added; we refer to Refs. [70] and [72] for a discussion on the coupling constants in the quark Lagrangian. Once the Lagrangian density is established, the variation of the action (integral of the Lagrangian density) with regards to the wave functions and fields in play, yields a set of equations that must be solved numerically.
Another phenomenological approach to compute the EoS are Skyrme density functionals. NS matter is out of the reach of perturbative Quantum Chromo-Dynamics (QCD), as is oftentimes represented in QCD phase diagrams. However, mesons that mediate the strong interaction can be treated in an effective-field theory designed to replace the extension of QCD to low temperature neutron matter. The Skyrme force is a nonrelativistic approach to the interaction of nucleons; the Skyrme density functional derived from this interaction is treated with the variational principle to design a Hamiltonian; baryons emerge as solutions to an approximation of Schrödinger equations (e.g., Hartree-Fock). In this framework, a series of (Skyrme) parameters (x i , with i ∈ [0, n]; t j , with j ∈ [1, m]; W 0 and α k with k ∈ [1, l]) are used. The upside of this phenomenological approach is that some microscopic quantities (energy per baryon, effective nucleon mass, symmetry energy, incompressibility etc.) are established analytically from the above-mentioned parameters; for details, see the review [78]. Nuclei data tables and properties of homogeneous neutron matter are used to adjust the parameters. Skyrme density functionals can be separated in two classes: standard and generalized. To classify Skyrme models, we refer to the expression for the Skyrme force ruling the interaction between nucleons presented in Ref. [79]: The quantity r ij is defined as the difference between the spatial coordinates of nucleon i and j, and p ij designates the relative momentum (difference between the momentum operator of i and of j). The spin-exchange operator between nucleons is denoted P s , and ρ(r ij ) is the local density or in other words the density at the barycenter r = 1/2(r i + r j ). The terms of Eq. (7) can be understood as follows: • terms proportional to t 0 are the effect of the force at zero range (hence the δ function), • terms proportional to t 1 and t 2 are effects for an effective range, and express the momentum dependence of the interaction -consequently finite temperature effects; let us note that t 4 and t 5 represent the effective range in the generalized form of the Skyrme force, and introduce density dependence to the term, • terms proportional to t 3 account for a three-body interaction expressed as a density-dependent two body interaction, • terms proportional to W 0 account for the two-body spin interaction with spin-orbit coupling.
To the Skyrme force can be added a pair force and Wigner force as is the case for Brussels-Skyrme models presented in this paper.
Although it is computationally costly, microscopic EoSs can be established by solving the N-body Schrödinger equation. In this framework, the only requirement to construct the EoS is a solid understanding of the nucleon-nucleon interaction (calibrated to nuclear data). Different approaches exist, for example, the nonrelativistic Brueckner-Hartree-Fock approach or the relativistic Dirac-Brueckner-Hartree-Fock [80] approach. In practice, the N-body problem is reduced to a three-body problem, which is sometimes itself reduced to a density dependent two-body problem. Another promising microscopic approach which has been largely explored for neutron-star matter in recent years is Chiral Effective Theory, but it only reaches densities that are twice the saturation density; in practice, one can glue low density Chiral Effective Theory results to a higher density EoS, as is presented in e.g., Ref. [81].
In this paper, we use the following cold matter EoSs: for references, see Table. I.
-DD2 and DDME2, -FSU2H, -H3, H4: based on the same nuclear model as the purely nucleonic GM1. The parameters for the nucleonic sector of this seminal model were chosen to reproduce simple nuclear constraints on the incompressibility (K = 300 MeV) and effective mass at saturation (m * /m = 0.7). Much more up-to-date and refined models are now available and used in this paper. We however employ GM1, H3, and H4 in this work as a comparison with [15]. Note that while GM1 and H4 have maximum masses larger than 2 M ⊙ , H3 does not: M max ≃ 1.79 M ⊙ , see Table. VI, -NL3, NL3ωρ.
• Hybrid relativistic mean-field (x5): for references, see Refs. [70] and [72]. -DD2-B15-40-20 -NL3ωρB20-50-0, NL3ωρ-B28-75-0 -NL3ωρ-B0-50-0  -NL3ωρ-B0-50-50 The terms for mesons self and cross interaction in the Lagrangian density are used to categorize the different family of models in Table I. Within a family of EoS, the difference between models lies in the value of the coupling constants. -DH: based on SLy4 parametrisation but contrary to the other SLy EoSs presented in this paper, it has been calculated consistently for the core and the crust; Classification of Skyrme EoSs and references are presented in Table II using the Skyrme parameters presented in Eq. (7).
• ab-initio (x1): BCPM as presented in Ref. [94]. It is based on the microscopic Brueckner-Hartree-Fock theory of nucleon interaction with the addition of a density dependent two-body force.
All above-mentioned EoSs with the exception of H3 meet the maximum mass criterion.

III. POLYTROPIC FITS
In the following, convention as presented in Ref. [15] (later referred to as PPFRead with PPF standing for Piecewise Polytropic Fits) is used: the rest-mass density ρ is directly connected to the baryonic density n via the baryon mass m B = 939 MeV/c 2 (ρ = n m B ). The energy density ǫ is calculated by using the first law of thermodynamics in the zero-temperature limit: A. Piecewise polytropes A polytrope takes the form with κ being the polytropic constant, and Γ being the adiabatic index. This type of crude approximation for the EoS is simple but oftentimes used: for example, the outer crust of NSs is well approximated by the pressure of ultrarelativistic electrons (Γ = 4/3). The whole NS cannot be accurately approximated with only one polytrope; however, it is possible to divide the EoS in N parts, each of which would pertain to a polytrope with fixed κ and Γ, which is what we call piecewise polytropes. A practical fit contains a restricted number of polytropes in order for the parametrisation to be convenient. PPFRead are based on N = 7 polytropes with four in the crust and three in the core and we shall use it as well. One could suppose that a total number of 3×N −1 parameters is needed to fit the EoS: N Γs, N κs and N − 1 transition densities ρ t[i→i+1] that defines areas of the EOS by which   (9) gives an expression for the energy density that depends on the adiabatic index and polytropic constant : with a i being the integration constant determined at the transition between polytropes. Expression for this constant can be found in Ref. [15], and its initial value is a physical requirement: the rest mass density zero limit implies ǫ = ρ. The major downside of piecewise polytropic fits is the nonderivability of quantities; for example, the sound velocity is not continuously defined.

B. Highlight on the importance of unified equations of state
It is common to find nonunified EoSs in the literature: the core and the crust EoS are not computed using the same nuclear model. However, nonunified constructions may induce non-negligible errors in the modeling of macroscopic parameters. Those constructions exist because computing the crust inhomogeneities is more tedious than computing the homogeneous core. A widespread practice within the astronuclear physics community, is to compute a core EoS, and glue an already established crust EoS. If this core-crust matching is performed with no care for thermodynamic consistency, it can lead to nonphysical jumps in the pressure. Only part of the outer crust is constrained by nuclei data, such that the high density outer crust and the inner crust are subject to model dependent differences. Therefore, core-crust matching of two EoSs with very different microscopic parameters such as the symmetry energy J m and its slope L m at the matching density, results in errors as large as 20% in the modeling of macroscopic parameters, for details see Ref. [13].
The relativistic hydrodynamics equations are nonlinear and the NS interior is opaque to observations. Therefore, one cannot separate the contribution of different parts of the EoS in the modeling of the total macroscopic parameters. Let us note that some quantities are more sensitive to certain parts of the star interior: for example the role of the crust treatment is more important for the radius of the star than for the mass, the total moment of inertia and the tidal deformability. If one wants to explore high density matter with NS observations, the artificial errors introduced by nonunified models can mislead into the acceptance or exclusion of the investigated nuclear model.
Examples for the use of nonunified EoSs are found in more simulations than can be listed: analytical representations used for GW data, finite temperature simulations, modelization of NS parameters in modified gravity, magneto-hydrodynamics, "universal" relations, etc.
All EoSs presented in Sec. II D are unified models: • The set of EoSs based on RMF models are taken from Refs. [69,77,107]; those that include hyperons are constructed consistently with the available experimental measurements of the properties of hypernuclei Refs. [69,76]. The EoS for the inner crust is calculated within the Thomas-Fermi approximation consistently with the EoS of the core [77,108]. The EoS of the outer crust has not been obtained consistently but taken from Ref. [109]. Several other outer crust EoSs are available (e.g., Refs. [110,111]) but they all are strongly constrained by nuclear physics data and therefore very similar. In addition, as mentioned in Ref. [107], we have checked that the use of an EOS for the outer crust not fully consistent with the rest of the star does not significantly affect the star properties, the radius in particular, for masses above 1.0 M ⊙ , which is the mass range of all NS currently observed.
• The set of EoSs based on Skyrme models are taken from Ref. [107]. The construction of the crust does not include shell effects and curvature terms which results in a mass shift with respect to experimentally measured masses, see Ref. [112]. Thus the EoS of the external part of the outer crust differs from the one we would get by employing experimental data but it is small enough to impact the relation between the mass and the radius by less than 1%.
• For the ab initio EoS, the outer crust is based on the formalism of Ref. [110]; a density functional designed from Brueckner-Hartree-Fock computations is used in the deformed Hartree-Fock-Bogoliubov formalism for nuclei not included in the data table [44]. For the inner crust, the energy density functional derived from Brueckner-Hartree-Fock calculations is used in the Thomas-Fermi approximation.
A few groups care for the unification of EoSs in a different manner, for example, Ref. [113] proposes a meta model paired to the extended Thomas Fermi approximation for the low density part of the EoS. References [114][115][116] use a set of causal model-independent EoS obtained from a Taylor expansion around the saturation density; at low density, a crust EoS is matched in a thermodynamically consistent way.
PPFRead are reduced to a core fit, to which the fit of EoS DH's crust is attached at a matching transition pressure which is not derived from either crust or core nuclear model. Their construction is the following: • the high density part of the EoS is fitted by three polytropes; • the low density part of the EoS is fitted by four polytropes based on EoS DH; • the point at which the high density and the lowdensity polytropes are matched is adapted to ensure a minimal fit error on the whole nonunified construction; • the last polytrope of DH is prolonged or shortened to ensure that it crosses the first polytrope of the high-density part.
Three models of our set of EoSs overlap with Read's fitted ones: DH, H3 and H4. We investigate the following constructions: • the unified EoS H3 and H4, • the unified EoS DH, • PPFRead for H3, H4 and DH, and present the relation between pressure and baryonic density in Fig. 6. Results are not presented for EoS H4 in this figure because our focus is on the lowdensity part of the star, and H3 and H4 diverge only in the core (the difference between H3 and H4 lies in the hyperonic meson couplings). The lowest-density parts of the EoSs, that is to say, ρ < 10 11.6 g/cm 3 , overlap for all constructions: this corresponds mostly to the outer crust which is calibrated to experimental data; therefore, it is similar for all nuclear models. Over 10 11.6 g/cm 3 , DH and unified H3 are different. The zoom in the figure highlights the matching area: for H3, the matching density is ρ H3 m = 7.9477 × 10 13 g/cm 3 ; for H4, ρ H4 m = 8.8774 × 10 13 g/cm 3 . In between the last point of convergence for all constructions, and the matching of DH and H3 PPFRead, the curves are different which attests to the limit of laboratory measurement calibrations of the low density part of the EoS. Although precautions are taken to avoid jumps in the pressure, the differences between the sole DH crust and the various core EoSs in PPFRead are not negligible.
From the above mentioned constructions, macroscopic parameters of a cold and non-rotating NS are computed using the relativistic equations of hydrostatics, the tidal deformability and the moment of inertia equations. Results for the radius R, dimensionless tidal deformability Λ, and moment of inertia I as a function of the mass M are presented in Fig. 7 for unified H3, H4, and DH as well as their respective PPFRead. The relative difference between the fits and unified tables for the variable X at given mass M is given by (11) and is presented in the right panels of Fig. 7. DH PPF-Read is constructed from DH low density polytropes and DH high density polytropes, therefore is unified: the red line in the right panel of Fig. 7 shows that it coincides with the tabulated unified EoS DH. We conclude that the fit method is powerful when applied to unified EoSs. In the case of PPFRead H3 and H4, the low density polytropes of DH are matched to the high density polytropes of H3 and H4 at ρ H3 m and ρ H4 m . Therefore, the whole fit is not unified and differs from our unified tables as is shown for green and blue lines in Fig. 7. This indicates that using nonunified EoSs for piecewise polytropic fits induces an artificial error on the macroscopic parameters.
From Fig. 7, the uncertainty related to the tidal deformability leads to two interesting points. On the one hand, despite the largest relative error being that of the tidal deformability, in the left panels of Fig. 7 the accuracy looks very similar for all three quantities. This is due to the fact that Λ is plotted on a logarithmic scale, because it spans over two orders of magnitude in the interesting range of masses. As a consequence, the relative change of Λ by ≅ 15% corresponds to ≅ 2% in mass and such error bars would be of a similar size in Fig. 7 (left  panel). On the other hand, one can notice that the sign of the relative difference in radius and that in tidal deformability are not the same. However, Λ is calculated according to the relation Λ = 2 3 k 2 C −5 [35,117], with k 2 being the tidal Love number solved simultaneously with the Tolman-Oppenheimer-Volkoff equations. The dimensionless deformability being proportional to k 2 R 5 , an increase in radius should correspond to an increase in the tidal deformability. The sign difference in the relative uncertainty of the two quantities in Fig. 7 can be explained by the large error of PPFRead with regards to the unified EoS on the quantity k 2 . It is particularly large for low mass stars, such that it dominates the R 5 factor. For higher mass stars, the uncertainty on k 2 is smaller and the scale of R 5 dominates. This large k 2 error can be understood as the strong effect of the crust matching on this quantity. Indeed, the fit from H3 unified tables in which the crust is treated correctly with the core -as is presented in details in the next section-gives a relative error on k 2 which is at most 3% (for a 1.5 M ⊙ ) whereas PPFRead give a relative error on k 2 of at most 50% (for a 1.0 M ⊙ ).
In the next section, we provide piecewise polytropic fits parameters based on unified tables of EoSs presented in Sec. II D.

C. Method for the fit
To revise PPFRead with unified tables of EoSs, we use an adaptive nonlinear least squares method to calculate the N poly − 1 transition densities ρ t , with the number of polytropes N poly = 7. They are adapted such that the fit error is minimized for the entire EoS (core and crust) of the P (ρ) fit.
For each unified table, we create a distribution for the number of density points ρ different from the original tabulated EoS. Our distribution of points allows us to give more importance to the accuracy of the fit in some parts of the star than others. We use a total of 1500 points, allocating 1/5 of the points in the crust and 4/5 of the points in the core, because it resulted in better fits of the full EoS. In each region, we distribute these points uniformly on a logarithmic scale. The fragment of the unified tables with the largest densities is not used because the density goes beyond the central density at maximum mass. We eliminate this highest density part by calculating the maximum mass of the star from the unified table, and only interpolate the EoS up to n max , thus increasing the fit accuracy for astrophysical quantities.
The core-crust transition density has a particular influence on the success of our fit method. This area is particularly sensitive to changes in the polytropic parameters and in turn, the points could be miss-allocated. Therefore, we chose to test values of the core-crust transition densities from n = 0.06 to n = 0.14 fm −3 for each EoS fit.
The relation P (ρ) is interpolated by using a first-order spline method to establish the pressure points from our distribution of density points. Each polytrope is fit by using a nonlinear least squares method from Eq. (9); Γ and κ are determined for each polytrope. Then, transition densities are recalculated from the polytropic parameters just fitted and the whole process with these new transition densities starts over until the set of ρ t stagnates. Finally, the energy density ǫ is calculated from Eq. (10).
Some alternatives to this fit method have been explored. With regards to the number of points for the interpolation, we have tested values between 200 and 10 000 points and observed a plateau of accuracy for ≅1500 points. We have also tested an inverse fit method, starting the fit from high density to low density which renders a similar accuracy.
The method described above is implemented to calculate the fit parameters which are presented in Tables IV, and V in the Appendix. We provide a routine in Python and tables of fit parameters in ASCII format to compute the tabulated fitted EoSs presented in this paper, see Ancillary files.

IV. ROLE OF UNIFIED EOSS IN THE ACCURACY OF MACROSCOPIC PARAMETERS
A. Piecewise polytropic fit's performance The method we use to establish piecewise polytropic fits of unified EoSs is intended to provide accurate modeling of macroscopic parameters for a nonrotating NS. Therefore, quantities M , R, I, and Λ are calculated from our fitted EoSs and compared with that of unified tables. Results for macroscopic parameters of prime interest are presented in Table VI, and VII: the maximum mass M max , density at the maximum mass n max , the radius at 1.0 M ⊙ and 1.4 M ⊙ respectively denoted R 1.0 and R 1.4 , as well as the radius at maximum mass R max , the moment of inertia at 1.338 M ⊙ denoted I 1.338 and at the maximum mass I max , and the tidal deformability at 1.4 M ⊙ denoted Λ 1.4 and at maximum mass Λ max . We provide the relative errors ∆ on these quantities except those defined at the maximum mass. For the latter we include instead the relative error δ defined as the relative difference between the quantities calculated at the maximum mass for the unified table, and at the maximum mass for our fit. Indeed since the maximum mass of the unified EoS and our fit are not equal, ∆ at the maximum mass and δ are different. The largest errors are presented in red in the tables.
With regards to nucleonic relativistic mean-field models, errors associated with the fit on M max , n max , and quantities related to the radius are systematically under under 1%. For astrophysical quantities related to the moment of inertia, it stays under 1.5%, and for the tidal deformability under 4%. For hyperonic relativistic mean-field models EoS consistent with the maximum mass criterion (all but H3), the errors associated to the fit on M max , n max , and quantities related to the radius stay under 1.5%. Quantities related to the moment of inertia stay under 2.5%, and for the tidal deformability under 7%. For hybrid relativistic mean-field models, the errors associated with the fit stays under 2% for all quantities.
For Skyrme models, the error associated with the fit for M max , n max , quantities related to the radius and the moment of inertia are under 1%. Once again, the tidal deformability does not fair as well, with an error up to 5%.
The maximum mass is the most accurately reproduced quantity, with an error under 0.5%. The tidal deformability is systematically the quantity with largest errors associated with the fit. Generally, nucleonic models are and H4 is presented in gray. For our fits, the relative difference is presented up to 98% of the maximum mass in plain lines, and the last two percent in dotted lines, see text for details. more accurately reproduced by our fits than hyperonic or hybrid ones. This is understandable because the number of polytropes chosen in the core is fixed to three, and the presence of hyperons or a phase transition to deconfined quarks in the core produces respectively an additional softening and drop of the adiabatic index as a function of the baryonic density (see Fig. 9 in the Appendix). In Fig. 8, the relative difference at given mass M between our fit and the unified tables for the radius, the tidal deformability and the moment of inertia are presented as a function of the mass for EoSs H3 and H4. The same quantities are plotted in gray for PPFRead to compare. The impact of using unified fits on the relative difference for the radius is particularly important for low mass stars, because the crust is relatively more significant for such objects. In overview, the relative difference for our fit is significantly smaller than for PPFRead for all macroscopic parameters considered.
The large increase of the relative error close to the maximum mass (Fig. 7 right panel and Fig. 8) is a conse-quence of choosing the mass as an independent variable. Close to the maximum mass point, the error increases significantly even for a small difference between the values of the maximum mass for the original and fitted EoS. This effect is absent for the dependence of the error as a function of central density or pressure. The relative inaccuracy of our fits for stellar configurations with the same central pressure is smaller in particular in the region close to maximum mass and for R, Λ and I are 0.1%, 1.8%, 1% respectively, compared with 1%, 8.7%, and 2% for a fixed mass M = 2 M ⊙ for EoS H4.
In overall, our fits perform well beyond the expected current precision from NICER for the radius and for the mass.

B. Universal relations
The story of "universal" relations started when [118] established the famous I-Love-Q (moment of inertia, tidal deformability and quadruple moment) relations; the quadruple moment Q expresses the deformation of the external gravitational field of a star by its rapid rotation (hence nonsphericity). It appeared that relations between some macroscopic parameters was only minimally EoS dependent. These relations have powerful predictive power, for the measurement of one parameter would permit the extraction of the other: when the first I-Love-Q relations were established, there was hope to extract the tidal deformability from the moment of inertia, when actually it is the other way around thanks to the detection of GWs. Many other relations appeared throughout the years and were used in various NS simulations. The physical meaning of this universality has been attributed to two main reasons: the first being the low density dependence of some of those relations and the second being the extension of the no-hair theorem 1 . The former refers to the calibration to laboratory experiments of the lowdensity part of the crust, but we have discussed that this argument is viable for part of the outer crust only. The latter argument is discussed in length in Refs. [120,121]: the authors explore the approximate baldness of NSs, and the analytical meaning of an extrapolated no-hair theorem in general relativity. It is suggested that this universality may rise from an emergent symmetry, acquired when the EoS parameters of a star are tuned out from main sequence, to relativistic stars, to black holes for which universality perfectly holds as per the no-hair theorem.
Universal relations have been established by fitting the modelled macroscopic parameters calculated from existing EoSs which are nonunified. It has been shown in Ref. [13] that the precision presented for some of those relations no longer held when compared with the modeled relation of parameters calculated with unified EoSs. Therefore, although the quasi-universality of those relations is not put into question, the precision of the fits can be tainted by the use of nonunified EoSs.
Recently, Godzieba et al. [122] proposed a revision of universal relations for multipole Love numbers based on the fit of polytropic type constructions. Under approximately the saturation density, the EoS is that of DH approximated by one polytrope, whereas, in the core, three polytropes are used with parameters adjusted to a Monte-Carlo-Markov-Chain: around two million EoSs are created, following basic rules of causality, maximum mass constraint around 2 M ⊙ , and GW170817 measurement of the tidal deformability; this construction of the EoS is not unified, however, authors use a method that adapts the matching density between [0.15-1.2] n 0 at the junction of DH and their polytropic cores, similarly to PPFRead. Universal relations are established for dimensionless electric tidal deformability Λ of order two, three and four with the compactness, but we focus on the relation C − Λ 2 which relates the compactness and the tidal deformability as measurable by GW detectors as of today. The fit is a logarithmic expansion as presented by [31] (Yagi & Yunes) and by Ref. [123] (Maselli et al.): The following EoS constructions are used: • unified EoS NL3 (stiff) and DD2 (soft), and BSR6 (in between), • DH crust matched to NL3, DD2, BSR6, H3 and H4 core at 0.15n 0 , • DH crust matched to NL3, DD2, BSR6, H3 and H4 core at 0.1 fm −3 , density at which [13] have shown that the uncertainty due to the use of nonunified model is minimized for models NL3, DD2 and BSR6, • DH crust matched to NL3, DD2, BSR6, H3 and H4 core at 1.2n 0 .
The  Fig. 10 for NL3 and DD2, and in Fig. 11 for H3 and H4. For both EoS NL3 and DD2, the quasi-universality of the relations is evident from the bottom plots of Fig. 10; all three universal relation fits overlap the relation C(Λ) for unified EoSs. In the case of H3 and H4 whose results are presented in Fig. 11 [13]. We find that the Yagi & Yunes fit is coherent with its 6.5% reported error, for all constructions of EoSs. In Ref. [122], the authors emphasize that their fit is an improvement on Yagi & Yunes one based on their collection of two million polytropic EoS constructions, however, when compared with tables of unified EoSs, the fit of Godzieba et al., which is based on nonunified constructions, generally fairs worse than the fit of Yagi & Yunes.
In overview, Yagi & Yunes fit performs the best out of the three fits. This fit falls within its reported precision when compared with unified as well as nonunified EoSs, for soft and stiff EoSs and modern or older (H3 and H4) EoSs.

V. CONCLUSION
We have established piecewise polytropic fits based on unified tables of equations of state for 52 different nuclear models: 15 nucleonic relativistic mean-field models, seven hyperonic relativistic mean-field models, five relativistic mean-field hybrid models, 24 nucleonic Skyrme models, and one ab-initio model.
The set of 52 equations of state chosen in this paper is confronted with both astrophysical observations and nuclear matter laboratory constraints. Microscopic constraints on nuclear models are discussed, particularly how the quantities J and L, respectively the symmetry energy and its slope calculated at saturation density, are constrained by laboratory measurements. The Direct Urca process, which is a very efficient cooling reaction in the core, occurs for 47 models of our set of equations of state.
We show that establishing piecewise polytropic fits based on nonunified equations of state results in errors on the modeling of macroscopic parameters that can be as high as approximately 15% ; note that this value corresponds to errors on the tidal deformability and is given excluding configurations close to maximum mass, because the relative uncertainty at fixed mass is then artificially high, as discussed in Sec. IV A. This artificial uncertainty is brought forth by the incompatibility of the crust and the core with regards to the microscopic quantities L and J.
We present parameters for piecewise polytropic fits based on unified equations of state. 14 parameters per model are required to establish the relation between the pressure, the energy density and the baryonic density; a Python routine and tables of fit parameters in ASCII format to compute the tabulated equation of state is provided, see Ancillary files. To establish the fit, a nonlinear least squares method is used, with adjusted transition densities. The precision of our fit was evaluated on macroscopic quantities.
For our 51 models consistent with M max 2M ⊙ , the fit error on key values of the mass, the radius, and the moment of inertia stays under 0.6% and 3.5% in the case of the tidal deformability for a broad range of NS masses except the region very close to the maximum mass, where this inaccuracy can be few times larger.
Finally, we confront the universal relations established by Ref. [122] between the compactness and the secondorder tidal deformability, to that of unified equations of state, as well as two other universal relations [31,123]. We conclude that the unified treatment of the crust plays a role in the reported precision of some of those relations. We present the parameters of the piecewise polytropic fits based on unified tables of relativistic mean-field models in Table IV and of Skyrme and ab initio models in Table V.

Appendix B: Accuracy on macroscopic parameters
We present the relative difference for astrophysical quantities of interest between the unified tables and our piecewise polytropic fits for relativistic mean-field models in Table VI and for Skyrme and ab initio models in Table VII.  V: Parameters of unified fits by seven polytropes of 24 nucleonic Skyrme and one ab-initio equations of state. The logarithm of the transition densities ρi between the polytropes is given in g/cm 3 . For each polytrope i, the adiabatic index Γi is presented. Only the first constant κ0 is presented because all others can be calculated from pressure continuity.
EOS log 10 (κ0) Γ0 log 10 (ρ1) Γ1 log 10 (ρ2) Γ2 log 10 (ρ3) Γ3 log 10 (ρ4) Γ4 log 10 (ρ5) Γ5 log 10 (ρ6) Γ6 (Nucleonic) Skyrme EOS BSk20 12.4732 1.6396 6.9219   in pourcent related to unified piecewise polytropic fit. The maximum mass (in solar mass) Mmax, the density (in fm −3 ) at maximum mass nmax, the radius (in km) for a M⊙ NS R1.0, the radius for a 1.4 M⊙ NS R1.4, the radius at maximum mass RM max , the moment of inertia (in 10 45 g.cm 2 ) for a 1.338 M⊙ NS I1.338 as measured in the double pulsar PSR J0737−3039, the moment of inertia at maximum mass IM max , the tidal deformability for a 1.4 M⊙ NS Λ1.4, and the tidal deformability at maximum mass ΛM max are presented. In red, we indicate the equations of state that gives the largest relative fit error in each category.