The nuclear symmetry energy plays a major role on the structure of neutron-rich finite nuclei as well as on the bulk properties of neutron stars (Li et al., 2014; Baldo and Burgio, 2016). In the past decades, several laboratory experiments (Russotto et al., 2023) have been performed in order to investigate the symmetry energy in finite nuclei, e.g., measurements of the nuclear masses (Möller et al., 2012), the nuclear dipole polarizability (Roca-Maza et al., 2015), the giant and pygmy dipole resonance energies (Klimkiewicz et al., 2007; Carbone et al., 2010), isobaric analog states (Danielewicz and Lee, 2014), and the neutron skin thickness (Adhikari et al., 2021; 2022). Several neutron stars (NS) properties are sensitive to the symmetry energy, e.g., its internal composition, the crust-core transition density and therefore the crust thickness, and the presence of fast direct URCA neutrino processes which regulate NS cooling (Yakovlev and Pethick, 2004; Burgio et al., 2021).

The symmetry energy is directly related to the more general and comprehensive task of the study of the equation of state (EoS), which plays a major role in nuclear structure studies, analysis of the heavy-ion collision dynamics, and the physics of compact objects (Oertel et al., 2017; Burgio and Fantina, 2018). The central density likely reached in NS interiors is about one order of magnitude larger than the nuclear saturation density, and this poses serious problems in theoretical astrophysics, because a correct theory of nuclear interactions for highly dense matter should be derived from the quantum chromodynamics (QCD). The well-known sign problem of lattice QCD still bars access to the high-density EoS, and therefore, models extracted from the nuclear many-body theory are required in order to build the EoS. Predictions have to be tested both in terrestrial laboratories, and in astrophysical observations. The most promising NS observables are the mass and radius; as far as the masses are regarding, the ones of several NSs are known with good precision (Lattimer, 2012; Demorest et al., 2010; Antoniadis et al., 2013; Fonseca et al., 2016; Cromartie et al., 2019; Romani et al., 2022), while the information on their radii (Özel and Freire, 2016; Guillot et al., 2013) has been improved thanks to the combined observations of NICER (Riley et al., 2019; Miller et al., 2019) and Advanced LIGO and Virgo collaborations, with the detection of gravitational waves emitted during the GW170817 NS merger event (Abbott et al., 2017; Abbott et al., 2018; Abbott et al., 2019). This event has provided us with important new information on the NS mass and radii by means of the measurement of the tidal deformability (Hartle, 1967; Flanagan and Hinderer, 2008), thus deducing upper and lower limits on it (Abbott et al., 2018; Radice et al., 2018). Further constraints on mass and radius have been recently reported by NICER for PSR J1231-1411, having mass M = 1.0 4 − 0.03 + 0.05 M ⊙ (Salmi et al., 2024).

In this paper, we concentrate on the study of the neutron skin thickness δ R in neutron-rich nuclei, such as ²⁰⁸Pb and ⁴⁸Ca, which has long been recognized as being strongly dependent on the slope of the symmetry parameter L . Novel data on ²⁰⁸Pb (PREX-I and PREX-II) (Adhikari et al., 2021) and ⁴⁸Ca (CREX) (Adhikari et al., 2022) with direct measurements consisting of parity-violating and elastic electron scattering technique (Horowitz et al., 2001), recently became available from the Thomas Jefferson Laboratory. Correlations between the neutron skin thickness, symmetry parameters, and NS observables like the radius of 1.4 M ⊙ and tidal polarizability have been widely analyzed, see, e.g., (Lattimer, 2023), and references therein.

In our previous paper (Burgio and Vidaña, 2020), we studied those kind of correlations by choosing a set of EoS based on microscopic methods and phenomenological approaches, and discussing their behaviour with respect to the PREX-I experimental data, which were the available ones at that time. Now, we would like to elaborate more on that analysis, taking into account the recent PREX-II (Adhikari et al., 2021) and CREX data (Adhikari et al., 2022). Moreover, we now choose a set of equations of state EoS which are compatible with the NS data on the highest observed mass M > 2.1 4 − 0.09 + 0.10 M ⊙ and tidal polarizability of a 1.4 M ⊙ , Λ 1.4 = 19 0 − 120 + 390 . This way, we should be able to improve our study on the possible correlations among observational quantities and properties of nuclear matter close to saturation.

The paper is organized as follows. In Section 2 we illustrate some basic properties of the EoS adopted in this work, along with the criteria selection for the choice of the optimal EoS. The laboratory and observational constraints on the nuclear EoS are presented in Section 3. The neutron skin thickness is discussed in Section 4, and conclusions are drawn in Section 5.

The composition of high density nuclear matter currently represents one of the most intriguing issues in theoretical physics, and several and diverse predictions have been proposed thus far (Burgio et al., 2021). The description of the extreme density conditions can include different scenarios, e.g., a purely nucleonic one characterized by a large neutron-proton asymmetry, or hyperonic matter or a hadron-quark phase transition. All those issues suffer of drawbacks that the current experimental data, either heavy-ion collisions in terrestrial laboratories or NS observations, cannot solve. In this work, we assume that nucleons are the most relevant degrees of freedom.

Theoretical approaches to determine the nuclear EoS are usually classified in microscopic and phenomenological ones. The interested reader is referred to recent reviews (Burgio and Fantina, 2018; Burgio et al., 2021); in this paper, we skip details and summarize the main features of the adopted methods. For the microscopic approaches, we adopt several EoS derived in the Brueckner–Hartree–Fock (BHF) theory (Baldo, 1999), which is based on the use of realistic two- and three-body forces (TBF), derived from meson-exchange theory (Machleidt et al., 1987; Nagels et al., 1978) and describe correctly the nucleonic phase shifts and the properties of the deuteron. For the TBF we use the phenomenological Urbana model (Pudliner et al., 1995; 1997; Baldo et al., 1997), and a microscopic TBF (Grangé et al., 1989; Baldo et al., 1997; Zuo et al., 2002; Li et al., 2008; Li and Schulze, 2008). We adopt as nucleon-nucleon potentials the Bonn B (Machleidt et al., 1987; Machleidt, 1989), the Nijmegen 93 (Nagels et al., 1978; Stoks et al., 1994), and the Argonne V 18 (Wiringa et al., 1995), which are supplemented by microscopic TBF and labeled in the following as BOB, N93 and V18. The Urbana model has been used with the Argonne V 18 potential and is labeled as UIX. The explicit inclusion of the quark-gluon degrees of freedom in the construction of a potential model has been performed in Ref. (Baldo and Fukukawa, 2014; Fukukawa et al., 2015), in which case two different EoS versions labeled respectively as FSS2CC and FSS2GC in Table 1 have been obtained. Besides BHF EoS, in this paper we exploit the often-used results of the relativistic Dirac-BHF method (DBHF) (Gross-Boelting et al., 1999), which employs the Bonn A potential, the variational APR EoS (Akmal et al., 1998) based on the V 18 potential, and the so-called CBF-EI model, obtained within the correlated basis function approach (Benhar and Lovato, 2017), using a realistic nuclear Hamiltonian with the Argonne V6’ (Wiringa and Pieper, 2002) and the Urbana IX nuclear potentials as TBF. For completeness, we also include in our set a parametrization of the Auxiliary Field Diffusion Monte Carlo (AFDMC) calculation (Gandolfi et al., 2010).

The philosophy of the phenomenological approaches is quite different from the one characterizing microscopic methods. In fact, they are based on effective interactions that are built to describe finite nuclei in their ground state, and therefore, predictions at high isospin asymmetries and density have to be taken with care (Stone and Reinhard, 2007). Among the most used ones, we mention Skyrme interactions (Vautherin and Brink, 1972) and relativistic mean-field (RMF) models (Boguta and Bodmer, 1977). In this work, we use a set of modern Skyrme EoS, which are listed in Table 1; in particular we mention the SLy0-SLy10 (Chabanat, 1995) and SLy230a (Chabanat et al., 1997; 1998) of the Lyon group, and the BSk20, BSk25 and BSk26 of the Brussels group (Potekhin et al., 2013; Goriely et al., 2013), the latter ones being unified EoS constructed on the basis of the energy-density functional theory. A complementary approach is given by RMF models, which are based on effective Lagrangian densities. The interaction between baryons is described in terms of meson exchanges, which are regulated by coupling constants of nucleons with mesons, and are usually fixed by fitting the bulk properties of nuclear matter as well as masses and radii of finite nuclei. In this work, we consider two types of RMF models: models with density-dependent coupling constants labeled DD-RMF (Nikšić et al., 2002; Typel and Wolter, 1999; Xia et al., 2022) and RMF models with fixed coupling strength (Mondal et al., 2016; Fattoyev et al., 2020; Das et al., 2021; Dhiman et al., 2007; Chen and Piekarewicz, 2015; Kumar et al., 2006; Sulaksono and Mart, 2006; Fattoyev and Piekarewicz, 2010; Kumar et al., 2017; Kumar et al., 2018; Routaray et al., 2023).

The main properties of the chosen EoS at saturation density are listed in Table 1. We notice that, whereas the saturation properties of the phenomenological models are within the empirical range, some microscopic EoS are marginally compatible with it. The reason is that the parameters of the phenomenological models are fitted on the saturation properties, while they are a prediction in the case of microscopic calculations, and those depend both on the many-body approach and the choice of the employed forces. For instance, the V18 EoS predicts a slightly too low saturation energy E 0 and compressibility K 0 , which is mainly due to the inclusion of a particular TBF (Li et al., 2008). We stress that a complete ab-initio theory of TBF is not available yet.

For completeness, we remind that the above mentioned methods are suited for describing the homogeneous component of the nuclear matter EoS, and that at densities ρ < ρ t ≈ 0.08 f m − 3 , clusters have to be included for the description of the NS crust. For that, we use the well-known Negele-Vautherin EoS (Negele and Vautherin, 1973) in the density range ( 0.001 f m − 3 < ρ < 0.08 f m − 3 ) , and the ones by Baym-Pethick-Sutherland (Baym et al., 1971) and Feynman-Metropolis-Teller (Feynman et al., 1949) for densities ρ < 0.001 f m − 3 typical of the outer crust.

An important check for the EoS is the behaviour of the symmetry energy slope L vs. S 0 , and this is plotted in Figure 1. The full green triangles represent the microscopic calculations (left panel), whereas the phenomenological ones are shown in the right panel. The experimental constraints indicate those derived from the study of isospin diffusion in heavy ion collisions (HIC, blue band) (Tsang et al., 2009); the electric dipole polarizability (violet band) (Roca-Maza et al., 2015); the neutron skin thickness in Sn isotopes (orange band) (Chen et al., 2010); the finite-range droplet mass model calculations (FRDM, magenta rectangle) (Möller et al., 2012); the isobaric analog state (IAS) phenomenology combined with the ²⁰⁸Pb neutron-skin thickness (green band) (Danielewicz and Lee, 2014); the recent analysis of the PREX-II experiment (black cross) (Essick et al., 2021). The blue solid curve is the unitary gas bound (Tews et al., 2017): only values of ( S 0 , L ) to the right of the curve are allowed. We see that all considered constraints are not simultaneously fulfilled in any area of the parameter space, probably because of the strong model dependencies in the extraction of the constraints from the raw data. Hence, at the moment, no theoretical models can be ruled out a priori, except those which are predicting values of the symmetry energy parameters outside the considered range.

A further important check regards the high-density behaviour of the nuclear symmetry energy, as illustrated in Ref. (Russotto et al., 2023). In the last few years several heavy-ion collisions experiments at relativistic energies have been performed in order to constrain the high-density symmetry energy. Figure 2 displays the ASY-EOS data (Russotto et al., 2016) (blue band) and the FOPI-LAND ones (Russotto et al., 2011) (light green band) as a function of the density, HIC (Sn + Sn) diffuseness measurements (Tsang et al., 2009) (grey band), whereas the red dashed contour labeled by IAS shows the results of Ref. (Danielewicz and Lee, 2014). For completeness, we also display the results of a Bayesian analysis (Tsang et al., 2024) which determines the boundaries at 68% (dark pink) and 95% confidence intervals (light pink) of the posterior distributions using an initial sample size of 3M of EoS. The experimental data are plotted up to ρ = 2 ρ 0 , and they all show a monotonically increasing behaviour with increasing density. In the four panels the symmetry energy is plotted vs. the nucleonic density for the microscopic models (upper left), for some of the Skyrme models (upper right), RMF models (lower left), and DD-RMF models (lower right) listed in Table 1. Except a couple of cases, i.e., BSk26 and SLy10, all EoS agree with the experimental data and the Bayesian analysis, thus confirming the need of more accurate experiments in order to disentangle the various theoretical approaches.

The strong correlation between the neutron skin thickness and the slope parameter L at normal nuclear saturation density was shown first by Brown and Typel (Brown, 2000; Typel and Brown, 2001), and confirmed later by other authors (Steiner et al., 2005; Centelles et al., 2009; Horowitz and Piekarewicz, 2001; Horowitz et al., 2001; Furnstahl, 2002). A measurement of the thickness allows to establish an empirical calibration point for the pressure of neutron star matter at subnuclear densities, and coupled with a NS radius measurement could determine the pressure at supranuclear densities. In fact, the neutron skin thickness and the NS size originate both from the pressure of neutron-rich matter, hence are sensitive to the same EoS. Therefore, the Typel-Brown correlation would be helpful in establishing the pressure–density relationship over a wide range of densities inside neutron stars.

The neutron skin thickness can be defined as the difference between the neutron ( R n ) and proton ( R p ) root-mean-square radii: δ R = ⟨ r n 2 ⟩ − ⟨ r p 2 ⟩ . Since the microscopic approaches discussed before are not suited for the description of finite nuclei, we prefer to use a different approach based on Ref. (Vidaña et al., 2009), in which an estimation of the neutron skin thickness of ²⁰⁸Pb and ¹³²Sn was made following the suggestion of Ref. (Steiner et al., 2005). In this case δ R is calculated to lowest order in the diffuseness corrections as δ R ∼ 3 5 t , being t the thickness of semi-infinite asymmetric nuclear matter t = δ c ρ 0 δ c 1 − δ c 2 E s 4 π r 0 2 ∫ 0 ρ 0 δ c ρ 1 / 2 S 0 / E sym ρ − 1 E SNM ρ − E 0 − 1 / 2 d ρ ∫ 0 ρ 0 δ c ρ 1 / 2 E SNM ρ − E 0 1 / 2 d ρ . (8)

In that expression, E s is the surface energy taken from the semi-empirical mass formula equal to 17.23 MeV, r 0 is obtained from the normalization condition ( 4 π r 0 3 / 3 ) ( 0.16 ) = 1 , and δ c is the isospin asymmetry in the center of the nucleus taken as δ c = δ / 2 according to Thomas-Fermi calculations. For consistency, we use this same method also for calculating the thickness in the phenomenological approaches. In Figure 3, we show the results of our calculations and compare them with experimental bands regarding CREX (left panel, magenta) and PREX I + II (right panel, cyan). Those experiments yield for PREX I + II a neutron skin thickness δ R ( P b 208 ) = 0.283 ± 0.071 f m , whereas the measurement of the neutron skin of ⁴⁸Ca with the same technique gives smaller values, i.e., δ R ( C a 48 ) = 0.121 ± 0.035 f m , thus pushing L towards larger or smaller values respectively. Data are shown as a function of the parameter L . The linear increase of δ R with L is not surprising because the neutron skin thickness in heavy nuclei is determined by the pressure difference between neutrons and protons, and this is proportional to the parameter L , that is, P ( ρ 0 , δ ) ≈ L ρ 0 δ 2 / 3 . We notice that the theoretical predictions show some correlation between δ R and L , as indicated by the linear fits (solid line) and by the value of the correlation coefficient, r = 0.75 , in both cases. This is slightly smaller than the previous result shown in ref. (Burgio and Vidaña, 2020), where the chosen EoS set was not filtered with respect to the NS observational data. It has to be noticed that, in both cases, the microscopic calculations alone seem to lie on a curve with a slope different than the one of the phenomenological calculations. However, due to the small number of points available for the microscopic approaches, no firm conclusion can be drawn. We also notice that the same EoS is unable to reproduce both the CREX and PREX I + II data; for instance, among the microscopic approaches, the CREX data set is well reproduced by the FSS2CC, V18, and CBF-EI models, whereas DBHF, BOB, and UIX fall in the range of PREX I + II. Those EoS differ not only by the many-body technique adopted, but also by the nucleon-nucleon interaction. A similar behavior can also be found for the phenomenological models such as Rs, SkT4, BSR8, and GL97 respectively. This might indicate a possible tension between the experimental data and the current understanding of the EoS. Though, it has to be stressed that the present calculations of the neutron skin thickness are based on the concept of semi-infinite nuclear matter (see Equation 8), which could be inappropriate. However, a similar result has already been found in other works, using suitable methods for finite nuclei (Lattimer, 2023). Further laboratory experiments on medium size nuclei, or a re-analysis of the current data could help to clarify this point.

Finally, in Figure 4, we display a correlation matrix among the saturation properties shown in Table 1, with Λ 1.4 , M m a x and the neutron skin thickness for ⁴⁸Ca and ²⁰⁸Pb. The matrix confirms the weak correlation of δ R with the L parameter already shown in Figure 3. A weaker degree of correlation is found between Λ 1.4 and L (r = 0.55), whereas no evident correlation between Λ 1.4 and S 0 and K 0 is found.

In this work, we have analyzed the predictions of microscopic and phenomenological EoS for the neutron skin thickness δ R of ⁴⁸Ca and ²⁰⁸Pb, and compared with the recent experimental data, CREX and PREX I + II. We have used an ensemble of different EoS models, that includes microscopic calculations based on the (Dirac)Brueckner–Hartree–Fock theory, the variational method, and Quantum Monte Carlo techniques, as well as several phenomenological Skyrme, and RMF models. The chosen EoSs are compatible with the constraints imposed by laboratory data on saturation properties of finite nuclei, and observational data by NICER and GW170817 regarding the NS mass and Λ 1.4 .

We have found a linear correlation between the neutron skin thickness δ R of ⁴⁸Ca and ²⁰⁸Pb and the L parameter, as already pointed out by several authors using nonrelativistic and relativistic phenomenological models. A weaker linear correlation with the tidal deformability Λ 1.4 is evident.

The most important result of our analysis is that the same EoS cannot reproduce at the same time the CREX and PREX I + II experimental data. Therefore, those measurements do not allow us to select the most compatible EoS among the ones considered in this work. Future NS observations, along with planned experiments in existing facilities or in next-generation radioactive ion beam laboratories, are fundamental to provide more stringent constraints on the nuclear EoS, thus finally improving our knowledge of the extreme density matter conditions.