The Equation of State for the Nucleonic and Hyperonic Core of Neutron Stars

We reexamine the equation of state for the nucleonic and hyperonic inner core of neutron stars that satisfies the 2$M_{\odot}$ observations as well as the recent determinations of stellar radii below 13 km, while fulfilling the saturation properties of nuclear matter and finite nuclei together with the constraints on the high-density nuclear pressure coming from heavy-ion collisions. The recent nucleonic FSU2R and hyperonic FSU2H models are fine tuned improving the density dependence of pure neutron matter at subsaturation densities. The corresponding nuclear matter properties at saturation, the symmetry energy and its slope turn out to be compatible with recent experimental and theoretical determinations. We obtain the mass, radius and composition of neutron stars for the two updated models and study the impact on these properties of the uncertainties in the hyperon-nucleon couplings estimated from hypernuclear data. We find that the onset of appearance of each hyperon strongly depends on the hyperon-nuclear uncertainties, whereas the maximum masses for neutron stars differ by at most 0.1 $M_{\odot}$, although a larger deviation should be expected tied to the lack of knowledge of the hyperon potentials at the high densities present in the center of $2 M_\odot$ stars. For easier use, we provide tables with the results from the FSU2R and FSU2H models for the equation of state and the neutron star mass-radius relation.


INTRODUCTION
The equation of state (EoS) of matter inside neutron stars has received a lot of attention over the last decades (Lattimer & Prakash, 2004Oertel et al., 2017). Besides black holes, neutron stars (usually observed as pulsars) are the most compact known objects in the universe. Their bulk features, such as mass and radius, strongly depend on the properties of matter in their interior and, hence, on the EoS.
With regards to mass determinations, the most precise measurements of masses are clustered around the Hulse-Taylor pulsar of 1.4M (Hulse & Taylor, 1975). However, accurate values of approximately 2M have been determined very recently. This is the case of the PSR J1614-2230 of M = 1.97±0.04M (Demorest et al., 2010) and the PSR J0348+0432 of M = 2.01 ± 0.04M (Antoniadis et al., 2013).
In view of these findings and future observations, it is opportune to analyze whether theoretical models for the EoS of dense matter can satisfy both the 2M max-imum mass constraint and radii below 13 km. Moreover, these models should fulfill the saturation properties 1 of nuclear matter and finite nuclei (or atomic nuclei). On the one hand, in order to obtain small neutron star radii, a softening of the pressure of neutron matter, and hence of the nuclear symmetry energy, around 1-2 times saturation density n 0 (n 0 ≈ 0.16 fm −3 ) is required (Lattimer & Prakash, 2007;Tsang et al., 2012;Ozel & Freire, 2016). On the other hand, the total pressure should be stiff enough in order to sustain 2M neutron stars. Very few models can reconcile simultaneously both constraints (small radius and large masses) and, at the same time, produce a precise description of finite nuclei (Jiang et al., 2015;Horowitz & Piekarewicz, 2001a,b;Chen & Piekarewicz, 2015a;Sharma et al., 2015).
Furthermore, as density increases inside neutron stars, the transition from nuclear to hyperonic matter would be favoured energetically (Ambartsumyan & Saakyan, 1960). Thus, the EoS softens as new degrees of freedom, hyperons, appear (Glendenning, 1982) leading to smaller neutron stars masses, below the 2M observations. This is known as the "hyperon puzzle", whose solution requires a new mechanism to stiffen the EoS: stiffer hyperon-nucleon and/or hyperon-hyperon interactions, repulsive three-body forces with hyperons, new hadronic degrees of freedom that push the onset of appearance of hyperons to higher densities or the phase transition to quark matter below the hyperon onset (see Ref. (Chatterjee & Vidana, 2016) and references herein).
In a recent paper (Tolos et al., 2017) we have obtained the EoS for the nucleonic and hyperonic inner core of neutron stars by reconciling the 2M mass observations with the recent analyses of radii below 13 km for neutron stars. Moreover, we have fulfilled the saturation properties of nuclear matter and finite nuclei (Tsang et al., 2012;Chen & Piekarewicz, 2014) as well as the recent constraints extracted from nuclear collective flow (Danielewicz et al., 2002) and kaon production (Fuchs et al., 2001;Lynch et al., 2009) in heavy-ion collisions (HICs). The study was performed in the relativistic mean-field (RMF) theory for describing both the nucleon and hyperon interactions and the EoS of the neutron star core. Two models were formulated, denoted as FSU2R (with nucleons) and FSU2H (with nucleons and hyperons), based on the nucleonic FSU2 model of (Chen & Piekarewicz, 2014).
In the present paper, we update the parameters of our two models in order to improve the behaviour of the EoS of pure neutron matter (PNM) at subsaturation densities by avoiding possible instabilities in the low-density region. We determine the properties at saturation of the modified interactions and we compare our results for the symmetry energy and the slope of the symmetry energy to recent experimental and theoretical determinations, while providing predictions for the neutron skin thickness of the 208 Pb and 48 Ca nuclei. Finally, we reinvestigate the mass-radius relationships for the two models, and estimate the impact on the neutron star masses, radii and composition of the uncertainties in the hyperon-nucleon couplings.
The paper is organized as follows. In Sec. 2 we present the RMF model for the determination of the EoS in betaequilibrated matter. In Sec. 3 we show the newly calibrated nucleonic FSU2R and hyperonic FSU2H models. Then, in Sec. 4 we display the results for the mass-radius relationship for neutron stars and in Sec. 5 we estimate the impact on the stellar properties of the uncertainties in the hyperon-nucleon couplings. We finally summarize our results in Sec. 6. Tables with numerical data of the EoSs are provided in the Appendix.

THEORETICAL FRAMEWORK
In the covariant field theory of hadronic matter, the baryons are treated as Dirac particles that interact through the exchange of mesons (Serot & Walecka, 1986). The formalism has been in wide use over the last four decades for describing the properties of the nuclear EoS and of finite nuclei in a relativistic quantum framework. A contemporary formulation of the Lagrangian density of the theory (Serot & Walecka, 1986, 1997Glendenning, 2000;Chen & Piekarewicz, 2014) may be written in terms of the contributions from the baryons (b), leptons (l=e, µ), and mesons (m = σ, ω, ρ, and φ) as where Ψ b and ψ l stand for the baryonic and leptonic Dirac fields, respectively. The mesonic and electromagnetic field strength tensors are The isospin operator is represented by the vector I b . The strong interaction coupling of a meson to a certain baryon is denoted by g (with N indicating nucleon) and the electromagnetic couplings by q, while the masses of the baryons, mesons, and leptons are denoted by m. The coupling constants of the above Lagrangian encode in an approximate way the complicated nuclear many-body dynamics. The g σN and g ωN couplings of the isoscalar σ and ω mesons to the nucleon determine the energy per particle and density of the nuclear matter saturation point, and, thus, are instrumental for the ground-state properties of finite nuclei. The g ρN coupling of the isovector ρ meson to the nucleon is key for the nuclear symmetry energy. Essentially, the symmetry energy measures the energy cost involved in changing all the protons into neutrons in nuclear matter (Li et al., 2014). Therefore, the g ρN coupling impacts on the properties of heavy neutron-rich nuclei and of neutron stars. The Lagrangian density (1), moreover, incorporates self-interactions of the meson fields. The σ-meson self-interactions, with the κ and λ couplings, were introduced by (Boguta & Bodmer, 1977) and allowed for the first quantitatively successful descriptions of nuclear matter and finite nuclei within the relativistic theory. These couplings soften the EoS at moderate densities and allow one to obtain a realistic compressibility of nuclear matter (Boguta & Bodmer, 1977;Boguta & Stoecker, 1983) in agreement with the values extracted from experiments on nuclear giant resonances and heavy ion collisions. 2 The quartic self-coupling ζ of the vector ω meson was introduced by (Bodmer, 1991). The ζ coupling must be nonnegative to prevent abnormal solutions of the vector field equation of motion (Bodmer, 1991;Mueller & Serot, 1996). It then implies an attractive nonlinear interaction that softens the EoS for high densities (Bodmer, 1991), thereby directly affecting the structure and maximum mass of neutron stars (Mueller & Serot, 1996). Finally, a mixed interaction between the ω and ρ mesons, with the coupling Λ ω , modulates the density dependence of the nuclear symmetry energy-which is related to the pressure of neutron matter-and influences the neutron radius of heavy nuclei and the radii of neutron stars (Horowitz & Piekarewicz, 2001a,b).
The Dirac equations for the different baryons and leptons are obtained from the Lagrangian density (1) as where the quantities denote the effective masses of the baryons. Let us mention that only the time-like component of the vector fields and the third component of isospin have been written in Eq. (2) due to the assumption of rotational invariance and charge conservation. The field equations of motion of the mesons follow from the respective Euler-Lagrange equations, see for example (Serot & Walecka, 1986). Altogether, the theory leads to a set of coupled nonlinear field equations that involve strong couplings. The exact solution of these equations is extremely complicated if one attempts to quantize both the baryon fields and the meson fields. Physically, the baryons are the constituents of the nuclear medium, whereas the mesons are the carriers of the interaction between baryons. Thus, in order to be able to solve the equations of the theory, it is meaningful to replace the meson field operators by their expectation values, which then act as classical fields in which the baryons move. This approach is known as the relativistic mean-field theory (Serot & Walecka, 1986).
Denoting the meson mean fields in uniform matter as σ = σ ,ω = ω 0 ,ρ = ρ 0 3 , andφ = φ 0 , the mesonic equations of motion in the mean-field approximation for the uniform medium are where I 3b is the third component of the isospin of a given baryon, and we use the convention that for protons I 3p = +1/2. The quantities are, respectively, the scalar and vector densities for the b baryon. In terms of the baryonic and leptonic Fermi momenta, k F b and k F l , and of the respective Fermi energies the scalar and vector densities for the baryons and the vector densities for the leptons are expressed as With the above ingredients, one can compute the energy density and the pressure of the system. The energy density is given by where the energy densities of baryons and leptons take the expressions We note that in obtaining Eq. (8) for the energy density, the equations of motion (4) were used to rewrite the contribution to ε of b (g ωbω n b + g ρbρ I 3b n b + g φbφ n b ). Finally, the pressure can be computed using the thermodynamic relation where the baryonic and leptonic chemical potentials are given by The cores of neutron stars harbor globally neutral matter that is in β-equilibrium. Therefore, the chemical potentials and the number densities of the different particles in a neutron star core are related by the conditions where b i and q i denote, respectively, the baryon number and the charge of the particle i. These relations, the Dirac equations (2) for the baryons and leptons, and the field equations (4) for the mesonic fields σ, ω, ρ and φ, are to be solved self-consistently for a given total baryon density n. Once the chemical potential and the density of each species have been obtained at the given n, one can determine the energy density and pressure of the neutron star matter for each density.

MODELS FOR THE EQUATION OF STATE
From the Lagrangian (1), in (Tolos et al., 2017) we formulated the models FSU2R (with nucleons) and FSU2H (with nucleons and hyperons), with a motivation for accomodating massive enough stars and the new astrophysical measurements of small stellar radii within a self-consistent microscopic theory of the EoS for the core of neutron stars. Note that this type of approach is different from-but complementary to-the methods where the astrophysical and nuclear observables are mapped onto the EoS through piecewise parametrizations of the EoS (Raithel et al., 2016;Lattimer & Prakash, 2016;Ozel & Freire, 2016). To build our models we started from the nucleonic FSU2 model of (Chen & Piekarewicz, 2014) that reproduces heavy neutron star masses but was not constrained to radii. The condition of small stellar radii imposed a soft nuclear symmetry energy in the theory. We showed that the resulting FSU2R and FSU2H models, besides the mentioned astrophysical constraints, can successfully describe the properties of finite nuclei and conform to the constraints on the nuclear EoS from kaon production and collective flow in HICs (Fuchs et al., 2001;Lynch et al., 2009;Danielewicz et al., 2002).
In the present work we start by introducing a modification of the parameters of our models FSU2R and FSU2H of (Tolos et al., 2017) in order to refine the behavior of the EoS of pure neutron matter (PNM) in the region of subsaturation densities. We report the new version of the parameters in Table 1 (it should be mentioned that the form of the equations of motion remains the same irrespective of the specific values of the coupling constants). We have changed the value of the quartic isovector-vector coupling Λ ω of FSU2R and FSU2H from 0.05 in (Tolos et al., 2017) to 0.045. This has been done because it results in a symmetry energy that is a little stiffer than before and avoids a previous instability in the EoS of PNM for low subsaturation densities. As Λ ω has been changed, we have refitted accordingly the value of the coupling g 2 ρN between the ρ-meson and the nucleons to obtain the same good reproduction of binding energies and charge radii of finite nuclei as in (Tolos et al., 2017). The values of the other parameters of FSU2R and FSU2H are the same of (Tolos et al., 2017). Owing to the fact the Λ ω and g 2 ρN couplings only contribute in neutron-rich matter, the EoS of symmetric nuclear matter (SNM), composed of the same number of protons and neutrons, is identical to that of our models FSU2R and FSU2H in (Tolos et al., 2017).
We collect in Table 2 a few characteristic isoscalar and isovector properties at the nuclear matter saturation density n 0 for the present version of our models. In the work (Fortin et al., 2015), the authors derived the constraint 1.7 P (n 0 ) 2.8 MeV fm −3 for the pressure of neutron star matter at saturation density.   Table 2 Properties at saturation of the models FSU2R and FSU2H of this work. We show the saturation density (n0), energy per particle (E/A), compressibility (K), and effective nucleon mass (m * N /mN ) in symmetric nuclear matter, as well as the symmetry energy (Esym), slope of the symmetry energy (L), curvature of the symmetry energy (Ksym), and pressure of pure neutron matter (P PNM ) at n0. They deduced this constraint from the results of the microscopic calculations of PNM performed by (Hebeler et al., 2013) using chiral two-nucleon and three-nucleon interactions, which are in good agreement with the results by (Gandolfi et al., 2012) from Quantum Monte Carlo calculations with the Argonne v 18 nucleon-nucleon potential plus three-nucleon forces. A narrower range 2.3 P (n 0 ) 2.6 MeV fm −3 was estimated more recently by (Hagen et al., 2015) from ab initio calculations of nuclear systems with chiral interactions. While our models FSU2R and FSU2H of (Tolos et al., 2017), with PNM pressures at saturation of 2.27 and 2.06 MeV fm −3 , fulfill the constraint of (Fortin et al., 2015), they are somewhat below the constraint of (Hagen et al., 2015). Now, with the new parametrization of our models, we are able to obtain PNM pressures at saturation density of 2.44 MeV fm −3 in FSU2R and of 2.30 MeV fm −3 in FSU2H (see Table 2) that are consistent with both the predictions from chiral forces derived by (Fortin et al., 2015) and (Hagen et al., 2015).
The EoS of PNM of the present FSU2R and FSU2H models differs from the results we showed in (Tolos et al., 2017) almost only for densities in the low-density region n n 0 , where the pressure of PNM is a little higher now. However, above saturation density, the pressures of our current parameters and those of (Tolos et al., 2017) are very similar. Consequently, compared with (Tolos et al., 2017), one may anticipate that the predictions for masses and radii of neutron stars will not be drastically affected.
The slope parameter L of the symmetry energy, i.e., , has become a standard reference in the literature for characterizing the stiffness of the change of the nuclear symmetry energy E sym (n) with density. In our original version of the FSU2R and FSU2H models shown in (Tolos et al., 2017), the value of the symmetry energy at saturation density Slope of the symmetry energy (L) versus symmetry energy (Esym(n 0 )) at the nuclear matter saturation density for the models FSU2R and FSU2H discussed in text. The shaded regions depict the determinations from (Li & Han, 2013;Lattimer & Lim, 2013;Roca-Maza et al., 2015;Hagen et al., 2015;Oertel et al., 2017;Birkhan et al., 2017).
was E sym (n 0 ) = 30.2 MeV in both models, while the slope parameter was L = 44.3 MeV in FSU2R and L = 41 MeV in FSU2H. In the updated version of FSU2R and FSU2H of the present work, these properties become E sym (n 0 ) = 30.7 MeV and L = 46.9 MeV in FSU2R and E sym (n 0 ) = 30.5 MeV and L = 44.5 MeV in FSU2H (see Table 2). These values suggest a relatively soft nuclear symmetry energy. We have plotted in Fig. 1 the ranges for E sym (n 0 ) and L that have been estimated in several recent works through the analysis of a variety of nuclear data from terrestrial experiments, astrophysical observations, and theoretical calculations (Li &  and FSU2H have an overlap with the majority of these ranges. We would like to remark that this is an a posteriori result, because the predicted values of E sym (n 0 ) and L are the consequence (Tolos et al., 2017) of having adjusted the FSU2R and FSU2H parameter sets to reproduce neutron star radii of about 13 km, without sacrificing maximum masses of 2M nor the description of binding energies and charge radii of atomic nuclei. Hence, we interpret the reasonable agreement of our results with the multiple constraints in Fig. 1 as hinting at the plausibility of the existence of neutron stars with relatively small radii. The neutron matter EoS is also strongly related with the neutron distribution in atomic nuclei. Models with softer symmetry energies produce a thinner neutron skin ∆r np (difference between the rms radii of the neutron and proton density distributions) in nuclei (Alex Brown, 2000;Horowitz & Piekarewicz, 2001a). Unfortunately, neutron densities and neutron radii are poorly known to date because the distribution of neutrons in a nucleus is hard to measure. Our present FSU2R and FSU2H models predict a neutron skin thickness of 0.15 fm in the neutron-rich nucleus 208 Pb. This prediction is compatible with the range 0.13 ∆r np 0.19 fm for 208 Pb extracted in (Roca-Maza et al., 2015) from measurements of the electric dipole polarizability of nuclei, the value ∆r np = 0.15 ± 0.03 fm determined from coherent pion photoproduction in 208 Pb at the MAMI facility (Tarbert et al., 2014), and the value ∆r np = 0.302±0.177 fm from parity violating electron scattering on 208 Pb performed at JLab (Abrahamyan et al., 2012;Horowitz et al., 2012). 3 In the case of the lighter nucleus 48 Ca, we find a neutron radius of 3.55 fm with FSU2R and of 3.57 fm with FSU2H, and a neutron skin of 0.166 fm with both models. The prediction is in good accord with the ranges 3.47-3.60 fm for the neutron radius and 0.12-0.15 fm for the neutron skin of 48 Ca obtained in (Hagen et al., 2015) through ab initio calculations of the neutron distribution of 48 Ca using nuclear interactions derived from chiral effective field theory; it also is in accord with the neutron skin of 0.14-0.20 fm for 48 Ca found from the new measurement of the electric dipole polarizability in 48 Ca (Birkhan et al., 2017). Altogether, it appears that the properties of the symmetry energy of the proposed models for the EoS, which are motivated by reproducing small neutron star radii (see next section), are compatible within uncertainties with different empirical and theoretical extractions of these properties. 3 We note that while experimental data are always provided with the associated error bars, theoretical models like ours, after the values of the coupling constants have been specified, make "exact" predictions with no error bars. In the future, it will be worth estimating error bars on our theoretical results, following recent initiatives to assess statistical errors and error propagation in nuclear functionals (Dobaczewski et al., 2014;Chen & Piekarewicz, 2014).

STELLAR PROPERTIES
Having access to the pressure and energy density of matter, we can compute the properties of neutron stars by solving the Tolman-Oppenheimer-Volkoff (TOV) equations (Oppenheimer & Volkoff, 1939). For static and spherically-symmetric stars, the TOV equations read as where r is the radial coordinate, m is the mass enclosed by a radius r, and G is the gravitational constant. For a given central density, the integration of these equations provides the corresponding mass and radius of the star. By repeating the calculation for different central densities, the mass-radius (M-R) relation of neutron stars can be obtained. Indeed, to solve the TOV equations for a neutron star we need the EoS of matter over a wide range of densities from the center to the surface of the star. The structure of a neutron star is such that the heavy liquid core is surrounded by a thin solid crust (Shapiro & Teukolsky, 1983;Haensel et al., 2006). The transition from the core to the crust occurs when the density of matter becomes lower than approximately 1.5 × 10 14 g/cm 3 . Below this density, matter ceases to exist in a homogeneous liquid phase because it is favorable that the protons concentrate with neutrons in nuclear clusters, which arrange themselves in a crystal lattice in order to minimize the Coulomb repulsion among them (Baym et al., 1971b,a;Shapiro & Teukolsky, 1983;Haensel et al., 2006). In the inner layers of the crust, the nuclear clusters are beyond the neutron drip point and the lattice is permeated by a gas of free neutrons in addition to the electron gas, whereas in the outer crust the nuclear clusters are neutron-rich nuclei below the neutron drip point, embedded in the electron gas. We have solved the TOV equations using the FSU2R and FSU2H models for the EoS of the uniform matter of the liquid core of the star for densities above 0.09 fm −3 (≈ 1.5 × 10 14 g/cm 3 ), under the conditions of β-equilibrium and global charge neutrality expressed in Eq. (12) of Sec. 2. At the densities of the crust, in the absence of calculations with our models of the complex structures that can populate this region of the star, we have used the EoS for the crust of neutron stars that has recently been derived from calculations based on the Brueckner theory in (Sharma et al., 2015).
The FSU2R model applies to nucleonic cores of neutron stars, i.e., when the whole stellar core is assumed to consist of neutrons, protons, electrons, and muons ('npeµ' matter). In the dense inner region of a neutron star core, however, the chemical potential may become so high that matter will be able to undergo a transition to other states of the low-lying octet of baryons, with hyperons appearing in the composition ('npYeµ' matter). Thus, we have devised the FSU2H model to allow for the presence of hyperons in the star interior (Tolos et al., 2017). In consequence, in the case of the FSU2H EoS, besides the nucleon and meson couplings shown in Table 1, we have considered the complete octet of baryons in the Lagrangian density (1). We have fixed the corresponding hyperon couplings from SU(3) flavor symmetry and from information on hyperon optical potentials in hypernuclei. We leave for the next section the discussion of the determination of the hyperon couplings and the analysis of the influence of the uncertainties associated with these couplings. Here, we focus on the results for the masses and sizes of neutron stars from the FSU2R nucleonic EoS and from the FSU2H hyperonic EoS with our baseline values for the hyperon couplings (given in Sec. 5).
We display the results for the relation between mass and radius of neutron stars in Fig. 2. A few data on the maximum mass configuration and the 1.5M configuration from FSU2R and FSU2H are presented in Table 3. For completeness, in Fig. 2, besides the curves of FSU2R and FSU2H, we also plot the M-R relations of two popular EoSs widely used in astrophysical calculations. They correspond to the Shen et al. EoS based on the relativistic TM1 nuclear mean field model (Shen et al., 1998) and to the Lattimer-Swesty EoS based on a non-relativistic Skyrme nuclear force (in its Ska version) (Lattimer & Swesty, 1991). Also shown is the result of the recent EoS from the Brueckner theory with the Argonne v 18 potential plus three-body forces computed with the Urbana model (Sharma et al., 2015). These additional EoSs are all non-hyperonic. We have included in the same Fig. 2 a few recent astrophysical determinations of neutron star mass-radius limits.
The M-R curve from each EoS exhibits a maximum mass, beyond which the star would become unstable against collapse into a black hole. The heaviest known masses of neutron stars are M = 1.97 ± 0.04M in the PSR J1614-2230 pulsar (Demorest et al., 2010) and M = 2.01 ± 0.04M in the PSR J0348+0432 pulsar (Antoniadis et al., 2013). We depict them by the horizontal bands in Fig. 2. Both the nucleonic FSU2R EoS and the hyperonic FSU2H EoS are able to provide maximum masses fulfilling the ≈ 2M observational limit, as well as the other EoSs shown in the same figure. We note that FSU2R reaches the maximum mass with a fairly compact stellar radius of 11.6 km (see Table 3). For canonical neutron stars with masses of 1.4M -1.5M , FSU2R predicts a radius of 12.8 km. The recent astrophysical determinations of neutron star radii from quiescent low-mass X-ray binaries in globular clusters and X-ray bursters seem to point in this direction  ., 2015)). The thin horizontal bands indicate the heaviest observed masses M = 1.97 ± 0.04M (Demorest et al., 2010) and M = 2.01 ± 0.04M (Antoniadis et al., 2013). The vertical blue band at the back depicts the M-R region constrained in (Hebeler et al., 2013) from chiral nuclear interactions up to n = 1.1n 0 and the conditions of Mmax > 1.97M and causality. The vertical red band at the front shows the M-R region derived from five quiescent low-mass X-ray binaries and five photospheric radius expansion X-ray bursters after a Bayesian analysis (Lattimer & Steiner, 2014b). The vertical striped yellow band is the M-R constraint derived from the cooling tails of type-I X-ray bursts in three low-mass X-ray binaries and a Bayesian analysis in (Nättilä et al., 2016) (model A of the paper). . Although these determinations are indirect and depend on stellar atmosphere models, they overall converge in favoring small neutron star radii in the range of about 9-13 km (Lattimer & Prakash, 2016;Ozel & Freire, 2016). An accurate radius measurement by new observatories such as NICER (Arzoumanian et al., 2014), which has begun operating aboard the International Space Station in June 2017, would represent a major step forward to corroborate or modify these expectations.
The compromise between having large maximum masses and small radii for canonical neutron stars is a challenging constraint that rules out a large number of theoretical EoSs (Lattimer & Prakash, 2016;Ozel & Freire, 2016;Oertel et al., 2017). This follows from the fact that the pressure of the high-density EoS must be hard enough to sustain massive stars, whereas the pressure at 1-2 times the nuclear saturation density n 0 must be, in contrast, effectively soft in order to pro- duce small radii for canonical mass stars. Given that the pressure of neutron star matter in the vicinity of n 0 is basically governed by the nuclear symmetry energy, the challenge is particularly acute in relativistic field theoretical models because the relativistic models usually have stiff symmetry energies. As we have demonstrated with FSU2R, it is possible to obtain parametrizations of the considered relativistic Lagrangian that meet both large stellar masses and-on condition of a soft symmetry energy-radii smaller than ∼ 13 km for M 1.4M , and that still provide an excellent reproduction of the binding energies and charge radii of finite nuclei (Tolos et al., 2017). There are some other parametrizations in the frame of the relativistic field theory that support these findings, such as the recent RMF012 and RMF016 models of (Chen & Piekarewicz, 2015b,a). Indeed, the accurately calibrated RMF016 model produces neutron stars of 2M and gives radii of 13 km for stars of 1.4M (Chen & Piekarewicz, 2015b,a), in keeping with the predictions of our FSU2R EoS.
When we allow for the appearance of hyperons in the neutron star core with the FSU2R model, the maximum mass of the star experiences a reduction of the order of 15%, due to the expected softening of the EoS, and then, with a maximum mass of 1.77M , it falls short of the 2M limit. In an effort to shed some light on the question whether with exotic degrees of freedom in the core, the star can satisfy the targets of 2M maximum mass and small radius at canonical mass, we have developed the hyperonic model FSU2H. In FSU2H, we essentially have stiffened further the nucleonic pressure above twice the saturation density, i.e., around the onset of appearance of hyperons. This comes at the price of a certain overpressure in symmetric nuclear matter for densities n 2n 0 when we compare it with the constraints deduced from the modeling of collective flow in HICs (Danielewicz et al., 2002), cf. Fig. 1 of (Tolos et al., 2017). Yet the pressure of FSU2H in pure neutron matter, shown also in Fig. 1 of (Tolos et al., 2017), fits within the projected region from the collective flow studies. Given that the β-equilibrated neutron-star matter is highly asymmetric, we consider this model as sufficiently realistic for describing neutron stars. The determination of narrower constraints on the EoS of PNM at several times n 0 from HIC experiments (Russotto et al., 2016) in the future should be of great help in this regard.
It can be observed in Fig. 2 that the FSU2H model with hyperons produces a comparable M-R relation to FSU2R and satisfies, as mentioned, the observational limit of 2M . With respect to FSU2R, in FSU2H the size of the radii has increased by 0.2-0.5 km for neutron stars heavier than 1M , expectedly, from the stiffer pressure of the nucleonic sector above twice the saturation density. The onset of hyperons occurs at a baryon density of 0.33 fm −3 , or 2.2n 0 . The maximum mass of 2.02M calculated with FSU2H is characterized by a radius of 12.1 km (see Table 3). For 1.5M stars, the hyperonic FSU2H EoS predicts radii of 13.2 km, which, although on the upper edge, are still compatible with the recent astrophysical indications of neutron star radii of about 9-13 km (Lattimer & Prakash, 2016;Ozel & Freire, 2016). The numerical results for the EoS and M-R relation of the FSU2R and FSU2H models are tabulated in the Appendix.
In closing this section, we ought to mention that the results for stellar radii of our EoSs have been possible while obtaining, within the same models, a realistic reproduction of the properties of atomic nuclei and of several other constraints. It seems unlikely that one may be able to account for significantly smaller neutron star radii in the theory considered here without abandoning the physical region of parameters. Hence, a discovery of even smaller stellar radii could provide evidence in favor of a phase transition to other degrees of freedom in neutron star interiors (Dexheimer et al., 2015).

IMPACT OF UNCERTAINTIES IN THE HYPERON COUPLINGS
We next discuss the determination of the values of the hyperon couplings in our FSU2H EoS and estimate the influence that the uncertainties in these couplings may have on the predictions for neutron star masses and radii.
The coupling of each hyperon to the scalar σ meson field is left as a free parameter to be adjusted to reproduce the hyperon potential in SNM, derived from hypernuclear data. It is well known that a Woods-Saxon type potential of depth U (N ) Λ (n 0 ) ∼ −28 MeV reproduces the bulk of Λ hypernuclei binding energies (Millener et al., 1988). As for the Σ hyperon, a moderate repulsive potential could be extracted from analyses of (π − , K + ) reactions off nuclei (Noumi et al., 2002) done in (Harada & Hirabayashi, 2006;Kohno et al., 2006). Fits to Σ − atomic data (Friedman & Gal, 2007) also point towards a transition from an attractive Σ-nucleus potential at the surface to a repulsive one inside the nucleus, the size of the repulsion not being well determined. The potential felt by a Ξ hyperon in SNM is also quite uncertain. Old emulsion data indicate sizable attractive values of around U (N ) Ξ (n 0 ) = −24 ± 4 MeV (Dover & Gal, 1983), while the analyses of the (K − , K + ) reaction on a 12 C target suggest a milder attraction (Fukuda et al., 1998;Khaustov et al., 2000). Taking these experimental uncertainties into account, we allow the hyperon potentials in SNM to take the following range of values: Note that we only consider uncertainties for the Σ and Ξ potentials, given the consensus on the Λ potential at saturation. The range of values for the hyperon potentials in SNM give rise to the following range for the hyperon-σ couplings: where the lower values correspond to the most repulsive situation (U The hyperon potentials are shown, as functions of the nuclear density, in Fig. 3. The left panel shows the potentials for isospin SNM, while the right panel corresponds to PNM, which is closer to the conditions of beta stable neutron star matter, where differences between the potentials for the different members of the same isospin multiplet can be seen. In order not to overcrowd the figure we have omitted the potentials of the positively charged hyperons, as they do not appear in the beta stable neutron star matter configurations appropriate for the present study. The coloured bands enclose the dispersion of results obtained employing the hyperonσ coupling ranges displayed in Eq. (17). The range of couplings has been determined from the uncertainties of the hypernuclear data and, strictly speaking, corresponds to normal nuclear matter density, n 0 . However, the coupling constants are density independent in our model and we can then obtain a range of values for the potential at any density. Specially important are the potentials around 2n 0 and beyond, which is the region of densities where hyperons are present in the models explored here (Tolos et al., 2017). As can be seen from the figure, the range of values for the hyperon potentials at 2n 0 in PNM are the following: U Ξ 0 (2n 0 ) = −22 to 5 MeV. We note that this range of values will strongly affect the composition of the neutron star, as we will show at the end of this section. In Fig. 4 we show the mass-radius relation for neutron stars obtained with the FSU2H model. The band collects the results obtained varying the hyperon couplings to the σ meson within the ranges in Eq. (17), which produce maximum masses that differ by at most 0.1M . This is a small effect, as it is obvious that the hyperon potentials at nuclear densities of around 6n 0 in the center of 2M stars (see Table 3) suffer a much larger uncertainty than the one we extrapolated from the normal nuclear densities characteristic of hypernuclear data. Indeed, the uncertainties tied to our lack of knowledge of the hyperon-nucleon and hyperon-hyperon interactions around the hyperon onset density of ∼ 2n 0 and beyond have often been exploited to build up RMF models that produce hyperonic neutron stars with maximum masses larger than 2M (Weissenborn et al., 2012;Bednarek et al., 2012;van Dalen et al., 2014;Oertel et al., 2015;Fortin et al., 2017). As can be seen in the extensive analyses of various models in (Fortin et al., 2015), the maximum masses turn out to be within a 0.3M band. It is therefore clear that determining the hyperon interactions at higher densities, as could be done from the analysis of HIC experiments (Morita et al., 2015), would help constraining the models in the appropriate regimes found in neutron stars.
Let us finish this section by showing, in Fig. 5, the particle fractions as functions of the baryonic density, for the FSU2H model (lower panel), where the coloured bands are obtained for the range of hyperon-σ couplings employed in this work. For completeness, we also show the particle fractions for the nucleonic FSU2R model in the upper panel, where we can see that the absence of negatively charged hyperons maintains a constant population of electrons and muons, and hence of protons and neutrons, already from slightly above 2n 0 . As for the FSU2H model, we note that all the particle fractions are affected by the hypernuclear data uncertainties, even if these are encoded only in the Σ and Ξ couplings to the σ meson. Upon inspecting the range of densities where hyperons may be present, we see that, although one can generally conclude that hyperons appear around 2n 0 , the order of appearance of each species is not determined, owing to the uncertainties derived from hypernuclear data. The first hyperon to appear can be either a Λ or a Σ − , the latter case only in the less repulsive situation allowed by data, namely when U (N ) Σ (n 0 ) 0 MeV. In fact, when the Σ feels its most repulsive potential value, it can even appear after the Ξ − hyperon. This happens when the Ξ potential value is on the most attractive side of the allowed region, namely U (N ) Σ (n 0 ) = −18 MeV. However, if one decreases the amount of attraction, as data permits, the Ξ − onset density is rapidly pushed towards larger values, even beyond the maximum density of 6n 0 represented in the figure, which stands as a representative central density of hyperonic stars. Summarizing, although hyperons are present in the interior of neutron stars modeled by the FSU2H interaction, the lack of precise knowledge on the hyperon-nuclear interactions prevents one from establishing the specific hyperonic composition in the interior of the star.
Using these updated interactions, we have obtained values for the PNM pressure at saturation density of 2.44 MeV fm −3 in FSU2R and of 2.30 MeV fm −3 in FSU2H, that are consistent with the estimates from chiral forces (Fortin et al., 2015;Hagen et al., 2015). The symmetry energy and its slope at saturation become E sym (n 0 ) = 30.7 MeV and L = 46.9 MeV in FSU2R and E sym (n 0 ) = 30.5 MeV and L = 44.5 MeV in FSU2H, thus being in good agreement with several recent estimates based on terrestrial experiments, different astrophysical observations, and theoretical calculations (Li & Han, 2013;Lattimer & Lim, 2013;Roca-Maza et al., 2015;Hagen et al., 2015;Oertel et al., 2017;Birkhan et al., 2017). Furthermore, the reviewed FSU2R and FSU2H models predict a neutron skin thickness of 0.15 fm in 208 Pb and of 0.166 fm in 48 Ca, which turn out to be compatible with previous experimental and theoretical determinations (Roca-Maza et al., 2015;Tarbert et al., 2014;Abrahamyan et al., 2012;Horowitz et al., 2012;Hagen et al., 2015;Birkhan et al., 2017) With regards to the mass and radius of neutron stars, radii below 13 km can be achieved because of the softening of the symmetry energy around saturation density whereas, at the same time, 2M stars can be obtained as the pressure of the high-density EoS is hard enough. These results are not drastically changed when using the updated FSU2R and FSU2H interactions as compared to the previous versions in Ref. (Tolos et al., 2017), because of the similar EoSs produced above saturation density. The numerical tabulations of the EoS and of the M-R relation from the models FSU2R (npeµ matter) and FSU2H (npYeµ matter) as a function of the number density n/n 0 are shown in Tables 4 and 5 for completeness.
However, the mass and composition of neutron stars might be strongly affected due to the uncertainties of the hyperon-nucleon couplings. The values of the hyperon couplings are determined from SU(3) flavor symmetry and from the available experimental information on hypernuclei, in particular by fitting to the optical potential of hyperons extracted from the data. The coupling of each hyperon to the σ meson field is left as a free parameter to be adjusted to reproduce the hyperon potential in SNM within the experimental uncertainties. As a result, we have found that the onset of appearance of the different hyperons strongly depends on the hyperon-nuclear uncertainties, whereas the maximum masses differ by at most 0.1 M , thus being less sensitive to the changes on the hyperon-nucleon couplings. This latter conclusion has to be taken with care, though, since the hyperon potentials at densities in the center of 2M stars suffer much larger uncertainties than the ones we have extrapolated from hypernuclear data at saturation. Hence, a greater dispersion of values for the maximum mass might be expected. The progress in the characterization of hyperon-nucleon interactions in dense matter derived from chiral effective forces (Haidenbauer et al., 2017), on the theoretical front, and from studies of HICs (Morita et al., 2015), on the experimental front, should contribute greatly to narrow down these uncertainties.

APPENDIX
For ease of use, in this appendix we provide in tabular form the results for the EoS and the M-R relation calculated with the FSU2R and FSU2H models discussed in the text. Table 4 Numerical data of the EoS for the core of neutron stars and of the M-R relation from the models FSU2R (npeµ matter) and FSU2H (npYeµ matter), as a function of the number density n/n0 (with n0 = 0.1505 fm −3 , cf. Table 2). The pressure P and mass-energy density ε are in MeV fm −3 , while the neutron star radius R and mass M are in km and M units, respectively.