A Modern View of the Equation of State in Nuclear and Neutron Star Matter

Abstract

Background: We analyze several constraints on the nuclear equation of state (EOS) currently available from neutron star (NS) observations and laboratory experiments and study the existence of possible correlations among properties of nuclear matter at saturation density with NS observables. Methods: We use a set of different models that include several phenomenological EOSs based on Skyrme and relativistic mean field models as well as microscopic calculations based on different many-body approaches, i.e., the (Dirac–)Brueckner–Hartree–Fock theories, Quantum Monte Carlo techniques, and the variational method. Results: We find that almost all the models considered are compatible with the laboratory constraints of the nuclear matter properties as well as with the largest NS mass observed up to now, $2.{14}_{-0.09}^{+0.10}{M}_{\odot }$ for the object PSR J0740+6620, and with the upper limit of the maximum mass of about 2.3–2.5$M_{\odot }$ deduced from the analysis of the GW170817 NS merger event. Conclusion: Our study shows that whereas no correlation exists between the tidal deformability and the value of the nuclear symmetry energy at saturation for any value of the NS mass, very weak correlations seem to exist with the derivative of the nuclear symmetry energy and with the nuclear incompressibility.

1. Introduction

A neutron star (NS) is the collapsed core of a massive star ($8-25M_{\odot }$, with $M_{\odot }\approx 2\times {10}^{33}$ g the mass of the sun), which at the end point of its evolution cannot be supported by hydrostatic pressure and collapses, producing a supernova explosion. A huge amount of gravitational energy is released, mainly in the form of neutrino radiation, and this leads to the complete destruction of the progenitor star. NSs may have masses in the range $M\sim 1-2M_{\odot }$ and radii of about $10-15$km. A huge amount of data has been collected from more than 50 years of NS observations, performed with ground-based and on-board telescopes covering all bands of the electromagnetic spectrum. In 2017, the observations witnessed an important breakthrough thanks to the direct detection (GW170817) of gravitational waves from such an event by the Advanced LIGO and Virgo collaborations Ref. [1,2,3]. In fact, it has been found that the observations of NS mergers can potentially provide strong constraints on the nuclear equation of state (EOS), as discussed in Refs. [4,5].

The EOS of isospin-asymmetric nuclear matter plays a major role not only in the study of NS structure and composition but also in the evolution of core-collapse supernovae and binary compact-star mergers [6,7]. Additionally, matter flows generated in heavy ion collisions (HIC) and the properties of nuclei in their ground state are strongly affected by the relevant features of the EOS, in particular the symmetry energy and its first derivative and the compressibility. In principle, it can be expected that in high-energy HICs as well as in supernova explosions and binary NS mergers, thermal effects are quite important, and thus they should be correctly included in the EOS. Besides that, the large density reached in the inner core of a NS can pose several theoretical problems because a complete theory of nuclear interactions based on QCD cannot be solved yet on the lattice for arbitrarily large values of density, temperature and isospin asymmetry, because of the well-known sign problem. Therefore, one has to rely on theoretical models and methods of the nuclear many-body theory in order to build the nuclear EOS.

Among possible observables regarding NSs, the mass and radius are the most promising ones, and they could ideally be used to infer the NS EOS within a certain observational uncertainty. While the masses of several NSs are known with good precision [8,9,10,11,12], information on their radii is less accurate [13,14]. However, simultaneous measurement of both quantities for several objects is required in order to constrain the EOS of NS matter and allow robust conclusions. The recent observation of gravitational waves (GWs) emitted during the merger of two corotating NSs [1,2,3] has opened the door to new possibilities of obtaining information on their masses and radii, by means of the measurement of the tidal deformability [15,16], and us allowed to deduce upper and lower limits on it [2,17]. The tidal deformability measures the linear response of the quadrupole deformation to a (weak) external gravitational field, and thus could be well constrained by the new data. It is therefore of interest to examine these quantities and their relations with other observables in theoretical calculations of the EOS. Moreover, the improved accuracy for the radius reached in recent observations by NICER (Neutron Star Interior Composition Explorer) Ref. [18,19], and the future planned missions like eXTP [20], will allow us to statistically infer NS mass and radius to better accuracy.

In this work we analyze the available constraints on the nuclear EOS, and compare with those derived from both phenomenological and ab-initio theoretical models. Possible correlations among the properties of nuclear matter close to saturation density with related quantities deduced from NS observations and nuclear physics experiments will be analyzed. We limit ourselves to the description of NS matter considering only nucleonic degrees of freedom, thus ignoring the possible appearance of hyperons [21] and a phase transition to the quark phase [22].

The article is organized as follows. In Section 2 we briefly review the different theoretical approaches for the nuclear EOS and illustrate the ones we adopt in the present study. We then discuss the experimental constraints on the nuclear EOS in Section 3. A brief overview of different EOSs of betastable matter is given in Section 4, and in Section 5 we discuss the comparison and correlations between NS EOS, nuclear physics and astrophysical constraints. Conclusions are drawn in Section 6.

2. Equations of State

The most commonly used theoretical approaches to determine the nuclear EOS can be classified into phenomenological and microscopic ones. (Non-)relativistic phenomenological approaches are based on effective interactions that are built to describe finite nuclei in their ground state, and therefore predictions at high isospin asymmetries should be considered with care [23]. In fact, at larger densities no experimental data are available, and therefore their behaviour can be very different. Skyrme interactions [24,25] and relativistic mean-field (RMF) models [26] are among the most used ones.

In this work we use a limited sample of Skyrme forces, namely GS and Rs [27], SLy4 [28] of the Lyon group, the old SV [29], SkI4 [30] of the SkI family, SkMP [31] and SkO [32]. We also include three EOSs derived within the modern Brussels-Montreal family of unified models, i.e., BSk22,24,26, which are commonly used in NS calculations [33]. We also consider two types of RMF models, which are based on effective Lagrangian densities where the interaction between baryons is described in terms of meson exchanges. In particular, we adopt models with constant meson-baryon couplings described by the Lagrangian density of the nonlinear Walecka model (NLWM), and models with density-dependent couplings [hereafter referred to as density-dependent models (DDM)]. Within the first type, we consider the models GM1 and GM3 [34]. For the DDM, we consider the models DDME1, DDME2 [35] and TW99 [36]. A further phenomenological RMF EOS, the SFHO EOS [37], has been used for comparison. A larger sample of RMF models, consistent with the analysis of Ref. [38], has been studied in Refs. [39,40], where the Love number and corresponding tidal deformabilities show very good agreement with the recent data from the GW170817 merger event.

Realistic two- and three-nucleon forces, that describe nucleon scattering data in free space and the properties of the deuteron, are the essential input for the microscopic approaches. These interactions are based on meson-exchange theory [41,42], or the recent chiral perturbation theory [43,44,45,46]. The main theoretical challenge is the treatment of the short-range repulsive core, which characterizes the nucleon-nucleon interaction, and this makes the difference among the available many-body approaches. The most well-known are the Brueckner–Hartree–Fock (BHF) [47] and its relativistic version, the Dirac–Brueckner–Hartree–Fock (DBHF) [48,49,50] theories, the variational method [51], the self-consistent Green’s function technique [52,53], the Quantum Monte Carlo techniques [54,55], the chiral effective field theory [56] and the $V_{\mathrm{low}}\mathrm{k}$ approach [57].

In this paper we adopt several BHF EOSs based on different nucleon-nucleon potentials, namely the Bonn B (BOB) [41,58], the Nijmegen 93 (N93) [42,59], and the Argonne $V_{18}$ (V18) [60]. In all those cases, the two-body forces are supplemented by nucleonic three-body forces (TBF), which are needed in all non-relativistic many-body methods in order to reproduce correctly the saturation properties of nuclear matter. Since a complete theory of TBF starting from first principles is not available yet, we adopt either phenomenological or microscopic models [61,62,63,64]. The phenomenological approach is based on the Urbana model (labelled as UIX) [62,65,66], whereas the microscopic TBF employes the same meson exchange as in the two-body force, as described in detail in Refs. [64,67]. Within the BHF framework, we also examine an EOS based on a potential model which includes explicitly the quark-gluon degrees of freedom, named FSS2 [68,69]. This reproduces correctly the saturation point of symmetric nuclear matter (SNM) and the binding energy of few-nucleon systems, and does not need TBF. We use two different EOS versions labelled respectively as FSS2CC and FSS2GC. Moreover, we compare these BHF EOSs with the often-used results of the Dirac-BHF method (DBHF) [49], which employs the Bonn A potential, in the following labelled DBHF(A), and a more recent calculation performed with the Bonn B potential [50], and labelled DBHF(B). We also compare with the APR EOS [51] based on the variational method and the $V_{18}$ potential, and a parametrization of a recent Auxiliary Field Diffusion Monte Carlo (AFDMC) calculation [70].

3. Bulk Properties of Nuclear Matter

Around saturation density ${\rho }_0$ and isospin asymmetry $\delta \equiv ({\rho }_n-{\rho }_p)/\rho =0$, being ${\rho }_n\left({\rho }_p\right)$ the neutron (proton) density and $\rho$ the total nucleonic density, the nuclear EOS can be characterized by a set of few isoscalar ($E_0,K_0$) and isovector ($S_0,L,K_{\mathrm{sym}}$) parameters, which can be constrained by nuclear experiments. The parameters are related to the coefficients of a Taylor expansion of the energy per particle of asymmetric nuclear matter as a function of density and isospin asymmetry,

$$\begin{array}{cc}\hfill E(\rho ,\delta )& =E_{\mathrm{SNM}}\left(\rho \right)+E_{\mathrm{sym}}\left(\rho \right){\delta }^2,\hfill \end{array}$$

(1)

$$\begin{array}{cc}\hfill E_{\mathrm{SNM}}\left(\rho \right)& =E_0+\frac{K_0}{2}{x}^2,\hfill \end{array}$$

(2)

$$\begin{array}{cc}\hfill E_{\mathrm{sym}}\left(\rho \right)& =S_0+Lx+\frac{K_{\mathrm{sym}}}{2}{x}^2,\hfill \end{array}$$

(3)

where $x\equiv (\rho -{\rho }_0)/3{\rho }_0$, $E_0$ is the energy per particle of symmetric nuclear matter at ${\rho }_0$, $K_0$ the incompressibility and $S_0\equiv E_{\mathrm{sym}}\left({\rho }_0\right)$ is the symmetry energy coefficient at saturation. These parameters are defined as

$$\begin{array}{cc}\hfill K_0& \equiv 9{\rho }_0^2\frac{d^2{E}_{\mathrm{SNM}}}{d{\rho }^2}\left({\rho }_0\right),\hfill \end{array}$$

(4)

$$\begin{array}{cc}\hfill S_0& \equiv \frac{1}{2}\frac{{\partial }^{2}E}{\partial {\delta }^2}({\rho }_0,0),\hfill \end{array}$$

(5)

$$\begin{array}{cc}\hfill L& \equiv 3{\rho }_0\frac{d{E}_{\mathrm{sym}}}{d\rho }\left({\rho }_0\right),\hfill \end{array}$$

(6)

$$\begin{array}{cc}\hfill K_{\mathrm{sym}}& \equiv 9{\rho }_0^2\frac{d^2{E}_{\mathrm{sym}}}{d{\rho }^2}\left({\rho }_0\right).\hfill \end{array}$$

(7)

The values of these parameters at ${\rho }_0$ for the various considered EOSs are listed in

Table 1. Saturation properties and NS observables predicted by the considered EOSs.

Model	EOS	${\mathit{\rho }}_0$[fm${}^{-3}$]	$-{\mathit{E}}_0$ [MeV]	${\mathit{S}}_0$ [MeV]	L [MeV]	${\mathit{K}}_0$ [MeV]	${\mathit{M}}_{\mathbf{max}}$[${\mathit{M}}_{\odot }$]	${\mathbf{\Lambda }}_{1.2}$	${\mathbf{\Lambda }}_{1.4}$	${\mathbf{\Lambda }}_{1.6}$
Micro.	BOB	0.170	15.4	33.6	70	238	2.50	1366	570	252
	V18	0.178	13.9	32.3	67	207	2.36	1082	442	188
	N93	0.185	16.1	36.5	77	229	2.25	1234	473	190
	UIX	0.171	14.9	33.5	61	171	1.96	848	309	112
	APR	0.159	15.9	33.4	51	233	2.20	720	274	110
	DBHF(A)	0.181	16.2	34.4	69	218	2.31	1635	681	295
	DBHF(B)	0.186	16.2	32.8	67	272	-	830	327	133
	FSS2CC	0.157	16.3	31.8	52	219	1.94	814	295	106
	FSS2GC	0.170	15.6	31.0	51	185	2.08	697	262	101
	AFDMC	0.160	16.0	31.3	60	239	2.21	822	293	109
Skyrme	Gs	0.158	15.6	31.2	94	239	2.13	1769	659	253
	Rs	0.158	15.1	30.8	86	248	2.12	1652	618	238
	SLy4	0.160	16.0	31.8	45	232	2.05	756	287	111
	SV	0.155	16.0	33.0	97	305	2.43	2224	914	393
	SkI4	0.158	16.2	33.7	106	245	2.17	1203	474	194
	SkMP	0.158	15.6	34.3	82	244	2.11	1295	487	188
	SkO	0.157	15.8	29.7	70	230	2.01	1252	451	164
	BSk22	0.158	16.1	32.0	69	246	2.26	1553	632	268
	BSk24	0.158	16.1	30.0	46	246	2.28	1260	523	227
	BSk26	0.159	16.1	30.0	38	241	2.17	830	327	133
NLWM	SFHO	0.157	16.2	32.8	53	244	2.06	862	334	132
	GM1	0.153	16.3	32.5	94	300	2.36	2223	913	393
	GM3	0.153	16.4	32.5	90	241	2.02	1688	617	228
DDM	DDME1	0.152	16.2	33.1	55	245	2.47	1765	773	355
	DDME2	0.152	16.1	32.3	51	251	2.51	1834	806	374
	TW99	0.153	16.2	32.8	55	240	2.08	1041	404	162
Exp.		∼ 0.14–0.17	∼ 15–17	28.5–34.9	30–87	220–260	$>2.{14}_{-0.09}^{+0.10}$	70–580
Ref.		[71]	[71]	[6,72]	[73,74]	[6,72]	[12]	[2]

.

From the measurements of nuclear masses [75] and density distributions [76] the values $E_0=-16\pm 1$ MeV and ${\rho }_0=0.14-0.17{\mathrm{fm}}^{-3}$ are obtained, whereas the value of $K_0$ can be extracted from the analysis of isoscalar giant monopole resonances in heavy nuclei. For the latter, results suggest $K_0=240\pm 10$MeV [77], or $K_0=248\pm 8$MeV [78], thus pointing to a rather soft EOS, as confirmed by HIC experiments [79].

Experimental information on the symmetry energy at saturation $S_0$ and its first derivative L can be obtained from the analysis of giant [80] and pygmy [81,82] dipole resonances, isospin diffusion measurements [83], isobaric analog states [84], measurements of the neutron skin thickness in heavy nuclei [85,86,87,88] and the meson production in HICs [89]. However, whereas $S_0$ is more or less well established ($\approx 3$MeV), the values of L ($30\mathrm{MeV}<L<87\mathrm{MeV}$), and especially those of $K_{\mathrm{sym}}$ ($-400\mathrm{MeV}<K_{\mathrm{sym}}<100\mathrm{MeV}$) are still quite uncertain and poorly constrained [90,91], and therefore we disregard them in our analysis.

From Table 1, we notice that all the adopted EOSs in this work agree fairly well with the empirical values, except the slightly too low $E_0$ and $K_0$ for V18, too large $S_0$ for N93, and too low $K_0$ for UIX and FSS2GC. We notice that several phenomenological models predict too large L values, whereas all the microscopic EOSs are largely compatible. This is clearly displayed in Figure 1, where L is plotted as a function of the symmetry energy at saturation $S_0$. The shaded areas represent the experimental data currently available. In particular, we report the constraints inferred from the study of isospin diffusion in HICs [92] (blue band), electric dipole polarizability [93] (violet band), the neutron-skin thickness in Sn isotopes [94] (grey region), the finite-range droplet mass (FRDM) model [95] and the isobaric-analog-state (IAS) phenomenology combined with the skin-width data (green diagonal region) [96]. Moreover the horizontal (red) band is obtained from a Bayesian analysis of mass and radius measurements of NSs [97], and the dashed curve is the unitary gas bound on symmetry energy parameters [90]. Only values of $(S_0,L)$ to the right of the curve are permitted, and therefore all the microscopic and some of the phenomenological models fulfill these constraints. We observe that there is no area of the parameter space where all constraints are simultaneously fulfilled, and this is likely due to the current uncertainties that plague the interpretation of the raw data. It should be stressed that the constraints on the EOS result from combining the raw data with theoretical models, and therefore they show some model dependence. On this basis, no theoretical model can be ruled out.

We report in Figure 2 the symmetry energy $E_{\mathrm{sym}}$ as a function of the baryon density. The results in the left (right) panel are plotted for the microscopic (phenomenological) EOSs, and are compared with the experimental data displayed by the shaded areas. In particular, the grey area represents the diffusion data of HICs, the green area includes the flow data obtained by the FOPI-LAND collaboration [98] on the collective flow, and the blue area is the experimental region checked by the ASY-EOS collaboration [99]. The full orange contour shows the results on the isobaric analog states (IAS), obtained in Ref. [96]. We see that most of the EOSs, both microscopic and phenomenological ones, are compatible with experimental data up to around the saturation density, whereas for larger densities some EOSs tend to predict smaller values for the symmetry energy that are below the experimental areas. This is a clear sign of discrepancy, which results in a much larger difference at larger values of the baryon density, such as the ones characterizing the inner core of a NS. We stress once again that the inferred constraints are model dependent, since the data interpretation requires theoretical simulations.

4. EOS for Betastable Matter

Once the EOSs for symmetric and pure neutron matter are defined, one can calculate the composition and the EOS of cold, neutrino-free, catalyzed matter. For charge-neutral matter in beta-equilibrium with neutrons, protons, and leptons ($e^{-}$, ${\mu }^{-}$), the EOS is computed in the following standard way [100]. One starts from the energy density of lepton/baryon matter as a function of the different partial densities ${\rho }_i$ of the species $i=n,p,e,\mu$,

$$\epsilon ({\rho }_n,{\rho }_p,{\rho }_e,{\rho }_{\mu })=({\rho }_n{m}_n+{\rho }_p{m}_p)+({\rho }_n+{\rho }_p)E({\rho }_n,{\rho }_p)+\epsilon \left({\rho }_e\right)+\epsilon \left({\rho }_{\mu }\right),$$

(8)

where $m_i$ are the corresponding masses, $E({\rho }_n,{\rho }_p)$ is the energy per particle of asymmetric nuclear matter, and the leptonic contribution is calculated assuming ultrarelativistic and relativistic expressions for the energy densities of electrons $\epsilon \left({\rho }_e\right)$ and muons $\epsilon \left({\rho }_{\mu }\right)$, respectively [100]. We have used the parabolic approximation [101,102] of the energy per particle of asymmetric nuclear matter in Equation (1), with the symmetry energy calculated simply as the difference between the energy per particle of pure neutron matter and symmetric nuclear matter,

$$E_{\mathrm{sym}}\left(\rho \right)\approx E({\rho }_n=\rho ,{\rho }_p=0)-E({\rho }_n=\rho /2,{\rho }_p=\rho /2).$$

(9)

From the energy density, Equation (8), the various chemical potentials can be computed,

$${\mu }_i=\frac{\partial \epsilon }{\partial {\rho }_i},$$

(10)

and imposing the beta-equilibrium conditions, ${\mu }_i=b_i{\mu }_n-q_i{\mu }_e$ ($b_i$ and $q_i$ denoting baryon number and charge of species i), along with the charge neutrality, ${\sum }_i{\rho }_i{q}_i=0$, one can find the equilibrium composition ${\rho }_i$ at fixed baryon density $\rho$, and finally the EOS,

$$p\left(\epsilon \right)={\rho }^2\frac{d}{d\rho }\frac{\epsilon \left({\rho }_i\left(\rho \right)\right)}{d\rho }=\rho \frac{d\epsilon }{d\rho }-\epsilon =\rho {\mu }_n-\epsilon .$$

(11)

Once the EOS of betastable matter is known, one can solve the Tolman-Oppenheimer-Volkoff (TOV) [100] equations which describe the structure of a non-rotating spherically symmetric star in general relativity

$$\begin{array}{cc}\hfill \frac{dp}{dr}& =-G\frac{\epsilon m}{r^2}\left(1+\frac{p}{\epsilon }\right)\left(1+\frac{4\pi p{r}^3}{m}\right){\left(1-\frac{2Gm}{r}\right)}^{-1},\hfill \\ \hfill \frac{dm}{dr}& =4\pi r^2\epsilon ,\hfill \end{array}$$

(12)

where G is the gravitational constant, p the pressure, $\epsilon$ the energy density and m the mass enclosed within a sphere of radius r. For each given central density, the integration of the TOV equations gives the mass and radius of the star corresponding to that density; this way one can construct an entire family of static configurations. It turns out that the NS mass has a maximum value as a function of the radius (or central density), above which the star is unstable against collapse to a black hole. The value of the maximum mass depends strongly on the nuclear EOS, hence the observation of a mass higher than the maximum mass allowed by a given EOS simply rules out that EOS.

We notice that the above mentioned theoretical methods cannot describe inhomogeneous and clusterized matter, and therefore for the low-density part $\rho <{\rho }_t\approx 0.08{\mathrm{fm}}^{-3}$, one has to adopt the well-known Negele-Vautherin EOS [103] in the medium-density regime ($0.001{\mathrm{fm}}^{-3}<\rho <{\rho }_t$), and the ones by Baym-Pethick-Sutherland [104] and Feynman-Metropolis-Teller [105] for lower densities $\rho <0.001{\mathrm{fm}}^{-3}$.

5. Constraints on the EOS from Terrestrial Laboratories and Astrophysical Observations

As already mentioned in Section 3, HICs at energies ranging from few tens to several hundreds MeV per nucleon have been exploited for extracting the gross properties of the nuclear EOS from the data. In fact, at a large enough energy, the two colliding nuclei produce flows of matter due to the large compression, resulting in a strong emission of nucleons and fragments of different sizes. The transverse flow, which is measured, depends sensitively on the pressure developed in the fireball at the moment of maximum compression during the collision. Additionally, the subthreshold $K^{+}$ production in heavy ion reactions has been demonstrated to probe the fireball density reached during the collision, with this being the ideal situation for exploring the EOS and its incompressibility.

In Ref. [106] the flow and kaon production analysis was summarized by plotting the region in the pressure versus density plane. A reasonable EOS should pass through it, and this is displayed in Figure 3 (left panels) as an orange area for the subthreshold kaon production [107], and as a grey area for the flow data [98]. Those results point in the direction of a soft EOS, i.e., values of the compressibility in the range $180\le K\le 250$MeV at density close to saturation. Those values are compatible with the ones extracted from the data on monopole oscillations [77]. In the upper (lower) panels of Figure 3 we plot the results for the microscopic (phenomenological) EOSs. We notice that most of the adopted EOSs are compatible with laboratory data, but some of them are too repulsive and therefore incompatible with experiments. In particular, among the considered microscopic EOSs, UIX, APR, and N93 are well compatible with the data extracted from HICs over the whole density range, whereas BOB, DBHF, and V18 are only marginally compatible large density (actually never reached in HICs), and those are characterized by a larger stiffness.

The EOS also rules the dynamics of NS mergers, in particular the final fate of the merger, a prompt or delayed collapse to a black hole or a single NS, as well as the amount of ejected matter which undergoes the nucleosynthesis of heavy elements, which strongly depends on the EOS. During the inspiral phase, the influence of the EOS is evident on the tidal polarizability, $\Lambda =\frac{2}{3}{k}_2{\beta }^{-5}$, where $k_2$ is the Love number and $\beta =GM/R$ is the compactness. An upper limit of $\Lambda <800$ was initially given in the first GW170817 analysis for a $1.4M_{\odot }$ NS [1], but later on the analysis was improved by assuming that both NSs have the same EOS, thus giving different limits of $\Lambda ={190}_{-120}^{+390}$, which translates into a measure of radius $R=11.9_{-1.4}^{+1.4}$km [2]. In this latter analysis, the values of the pressure as a function of density were extracted, and those are displayed as the blue hatched areas in Figure 3 (right panels). We notice that in this case the comparison has to be performed for the betastable case. We observe that almost all microscopic EOSs, except UIX, turn out to be compatible with the GW170817 data at density $\rho >2{\rho }_0$, whereas the nuclear collision data look more restrictive. Additionally, for the phenomenological case some EOSs turn out to be marginally compatible with the observational data, as for the flow data in the symmetric case.

A very important constraint to be fulfilled is the value of the maximum mass for the different EOSs, which has to be compatible with the observational data. In Figure 4 we display the mass-radius relations obtained with microscopic and phenomenological EOSs, shown respectively as solid and broken curves. We observe that most models give values for the maximum mass larger than $2M_{\odot }$, except the soft microscopic UIX and FSS2GC, which therefore are compatible with current observational data [9,10,11], in particular with the largest mass observed up to now, $2.{14}_{-0.09}^{+0.10}{M}_{\odot }$ at 68% confidence interval for the object PSR J0740 + 6620 [12] (dark orange band). For completeness, we also display the observational limits at 95% confidence interval (light orange band). Analogous limits are plotted also for the object PSR J0348 + 0432 [11] (grey bands), which are more restrictive at high confidence level. Apart from these lower limits, some recent theoretical analyses of the GW170817 event indicate an upper limit on the maximum mass of about $2.33M_{\odot }$ (68%) or $2.5M_{\odot }$ (95%) (displayed by red horizontal lines) [108,109,110,111], with which several of both the microscopic and phenomenological EOSs would be compatible as well. We also display the Bayesian parameter estimation of the mass and equatorial radius of the millisecond pulsar PSR J0030 + 0451 [18,19], as recently reported by the NICER mission. The $M,R$ values inferred from the analysis of the collected data (green and light grey zones) are $1.{36}_{-0.16}^{+0.15}{M}_{\odot }$ and $12.{71}_{-1.19}^{+1.14}$km [18], or $1.{44}_{-0.14}^{+0.15}{M}_{\odot }$ and $13.{02}_{-1.06}^{+1.24}$km [19].

Let us now turn to the discussion of the tidal deformability and its possible correlations with nuclear matter properties in its ground state. As already anticipated, the analysis of the GW170817 event [1,2,3] produced a value of $\tilde{\Lambda }<730$, assuming equal mass merging. If both NSs have the same EOS, this leads to the constraints $70<{\Lambda }_{1.4}<580$ and $10.5<R_{1.4}<13.3$km [2] for a $1.4M_{\odot }$ NS. A more stringent lower limit $\tilde{\Lambda }>400$ [17] on the average tidal deformability was imposed by the high luminosity of the kilonova AT2017gfo following the NS merger event. This constraint could indicate that $R_{1.4}\gtrsim 12$km [112,113,114,115], but it has to be taken with great care [116,117].

The possibility of finding correlations between properties of nuclear matter and NS observables has been recently explored [118,119]. In the following we further explore this issue, using the set of microscopic and phenomenological EOSs listed in Table 1. In Figure 5 we show the tidal deformability of a 1.2 (upper panels), 1.4 (central panels) and 1.6 (lower panels) solar-mass NS as a function of the symmetry energy at saturation $S_0$ (left panels), its first derivative L (central panels) and $K_0$ (right panels). The light- and dark-shaded bands in the central panels represent the limits inferred from the observational data of the GW170817 event [2] together with the experimental limits reported in Table 1. The degree of correlation is quantified by the correlation factor

$$r(x,y)=\frac{1}{n-1}\frac{{\sum }_x{\sum }_y(x-\overline{x})(y-\overline{y})}{{\sigma }_x{\sigma }_y},$$

(13)

with n being the number of data pairs, $\overline{x}$ and $\overline{y}$ being the mean values of x and y, and ${\sigma }_x$ and ${\sigma }_y$ being their standard deviations. We obtain the following values

$$\begin{array}{cc}\hfill r([S_0,L,K_0],{\Lambda }_{1.2})& =[0.006,0.635,0.709],\hfill \end{array}$$

(14)

$$\begin{array}{cc}\hfill r([S_0,L,K_0],{\Lambda }_{1.4})& =[0.206,0.551,0.702],\hfill \end{array}$$

(15)

$$\begin{array}{cc}\hfill r([S_0,L,K_0],{\Lambda }_{1.6})& =[0.145,0.459,0.682].\hfill \end{array}$$

(16)

As can be seen the correlation factor is rather small in the case of $S_0$ for the three values of the tidal deformability, indicating that no correlation at all exists between $\Lambda$ and $S_0$ for any value of the NS mass. Instead, a very weak correlation of $\Lambda$ seems to exist with L and a slightly stronger, although still weak, exists with $K_0$. We observe that in the case of ${\Lambda }_{1.4}$ all but two (SFHO and TW99) of the NLWM and DDME models are incompatible with the observational constraint [2]. On the contrary, most of the Skyrme models lie within the shaded bands, except for a few cases. Regarding the microscopic models, they are almost all in agreement with the GW observations and experimental constraints on $S_0$ and L, although only four are compatible with the constraints on $K_0$.

6. Conclusions

In this work, we have analyzed several constraints on the nuclear EOS currently available from NS observations and laboratory experiments. For this purpose, we have used a set of different models that include several phenomenological EOSs based on Skyrme and relativistic mean field models as well as microscopic calculations based on the (Dirac–) Brueckner–Hartree–Fock theories, the variational method and Quantum Monte Carlo techniques. To select the most compatible EOSs among the ones considered in this work, we have employed in particular the experimental constraints on several properties of nuclear matter at saturation density derived from different experiments as well as observational constraints on the mass, radius and tidal deformability imposed by recent measurements of the masses of millisecond pulsars [12], the data of the NICER mission [18,19] and the GW170817 NS merger event [1,2,3]. We have found that almost all considered models are compatible with the laboratory constraints of the nuclear-matter properties as well as with the largest masses observed up to now, $2.{14}_{-0.09}^{+0.10}{M}_{\odot }$ for the object PSR J0740 + 6620 [12], and with the upper limit of the maximum mass of about 2.3–2.5$M_{\odot }$ [108,109,110,111] deduced from the analysis of the GW170817 event. Our study of possible correlations among properties of nuclear matter at saturation density with NS observables, particularly with the tidal deformability, has shown that no correlation exists between $\Lambda$ and $S_0$ for any value of the NS mass, but weak correlations of $\Lambda$ do exist with L and with $K_0$.

We would like to finish by noticing that while the isoscalar part of the nuclear EOS is rather well constrained by the major experimental, observational and theoretical advances, the isovector one is less well known mainly due to our still limited knowledge of the nuclear force and, particularly, of its in-medium modifications and its spin and isospin dependence. Future NS observations, such as the precise simultaneous measurement of the mass and radius of a single NS, together with laboratory experiments planned in next-generation radioactive ion beam facilities, are fundamental to provide more stringent constraints on the nuclear EOS, and are very much awaited for.