The era of precision observation of the neutron star, and the multi-messenger astronomy-astrophysics has been opened. One of the most fundamental questions in the neutron star physics is what is the correct equation of state (EoS) of nuclear matter at densities below and above the saturation density (ρ₀). Equation of state is determined by the relation between the energy of a particle and its pressure which is gradient of energy with respect to density. Basic properties of the neutron star such as mass, size, density profile, tidal deformability, species of particles, and their fractions are determined directly by the EoS. Neutron star cooling is critically dependent on the fraction of particles, so the EoS also plays an essential role in the thermal evolution of neutron stars.

Energy of a particle in nuclear matter is conventionally expanded around the saturation density in terms of x = (ρ − ρ₀)/3ρ₀ as

Neutron-proton asymmetry is accounted by δ = (ρ_n − ρ_p)/ρ where ρ_n and ρ_p are the density of the neutron and the proton, respectively. Most of the nuclear structure models fit their model parameters to the nuclear data, and the constants and coefficients in Eqs.(2, 3) which determine the EoS of nuclear matter are obtained as results of the fitting. Among the constants and coefficients in Eqs.(2, 3), saturation density ρ₀ and the binding energy E_B are determined highly model-independently, and the values are converged to ρ₀ = 0.16 fm⁻³ and E_B = 16 MeV. On the other hand, compression modulus K₀ that characterizes the stiffness of EoS in the symmetric nuclear matter is uncertain much more than ρ₀ and E_B. For example, collection of 240 Skyrme force models shows that the models predict the K₀ value in the range 200 − 400 MeV [1]. Experimentally constrained range is reduced to 200 − 260 MeV, but it is still uncertain much more than ρ₀ and E_B. Since the density at the center of a canonical neutron star, a star of mass 1.4M_⊙ is in the range (2 − 3)ρ₀, role of K₀ in the neutron star matter might be non-negligible and therefore it is needed to reduce the uncertainty of K₀ in order to have more exact EoS of the neutron star matter. Situation is no better than K₀ for the symmetry energy S(ρ) in Eq. (3). Ranges of J, and L constrained from experiments are 30 − 35 MeV and 40 − 76 MeV, respectively [1]. For K_sym experimental access is still very limited, and the value is guessed from relations with L or 3J − L [2]. Properties of neutron rich systems (δ ≃ 1) are sensitive to the behavior of the symmetry energy. J, L and K_sym being siginificantly uncertain, large uncertainty in the EoS of neutron stars is unavoidable.

Goal of this work is to reduce the uncertainty of K₀, J and L so that we have more reliable EoS for the neutron star matter at densities up to (2−3)ρ₀ below which nucleons may be the dominant constituents of the matter. It has been shown that if the energy density of nuclear matter is expanded in powers of the Fermi momentum of the nucleon, 7 terms are sufficient for the correct description of EoS at densities relevant to nuclei and the center of neutron stars [3–5]. Higher order terms added to the basic 7 terms give vanishing contribution up to ρ ≃ 3ρ₀, and their effects are marginal at densities higher than 3ρ₀. Following this observation, we describe the EoS of infinite nuclear matter with 7 density-dependent terms.

Parameters of the KIDS functional are determined in two steps. In the first step, seven parameters in the nuclear matter functional are fixed to the values of ρ₀, E_B, K₀, J, L, K_sym, and Q_sym. After the 7 nuclear matter parameters are determined, the model is applied to nuclei. In the second step, energy per particle (E/A) and the charge radii (R_c) of spherical magic nuclei ⁴⁰Ca, ⁴⁸Ca and ²⁰⁸Pb are reproduced by adjusting two additional parameters that account for the gradient of density and spin-orbit interactions in nuclei [6]. At this stage the number of model parameters is increased to 9, but we don’t need to introduce additional parameters any more. Steps 1 and 2 are repeated by varying K₀, J, L and K_sym values (ρ₀, E_B, Q_sym kept unchanged), and χ² values are calculated for E/A and R_c of ⁴⁰Ca, ⁴⁸Ca and ²⁰⁸Pb for a given set of (K₀, J, L, K_sym). We find that nuclear data are most accurately reproduced in the range K₀ = 230−250 MeV. In the application to the neutron star, we select two sets of (J, L, K_sym) that give the smallest χ² values for a given K₀ value. Solving TOV equations, we obtain the mass and radius of the neutron star. Comparison with neutron star data provides us the ranges of K₀, J and L that are most consistent with both nuclear data and astronomical observation. As a result, we obtain K₀, J and L constrained in ranges K₀ ~ 230−250 MeV, J ~ 31−33 MeV and L ~ 55−65 MeV.

In Section II, we introduce the model and the strategy of investigation. Neutron star properties are considered, and reduced ranges of the EoS parameters are obtained in Section III. Parameters that are consistent with the neutron star observation are applied to calculating the properties of nuclei in Section IV. We summarize the work in Section V.

KIDS (Korea-IBS-Daegu-Sungkyunkwan) density functional theory provides nuclear energy density functional based on rules for systematic expansion of the energy of a nucleon in nuclear matter. In the homogeneous infinite matter, energy per particle is expanded in terms of the Fermi momentum as

where T is the kinetic energy, and the terms in the summation denote the potential energy originating from the strong interaction. In the cold nuclear matter kF ∝ ρ^1/3, so expansion in powers of ρ^1/3 is equivalent to the expansion in terms of the Fermi momentum. It’s been shown that 3 α_k’s and 4 β_k’s are suffcient for a correct description of the EoS of nuclear matter at densities well below and above the saturation density [5].

By transforming the KIDS functional to the Skyrmetype potential, one can obtain the single-particle potential of the nucleon in nuclei. In addition to α_k and β_k, two parameters accounting for the momentumdependent and spin-orbit interactions are introduced in the single-particle potential for the nucleon in nuclei. We have 9 free parameters altogether. Determination of the 9 parameters is divided into two steps: At first, we assume specific values of the 7 EoS parameters ρ₀, E_B, K₀, J, L, K_sym and Q_sym, and determine 3 α_k’s and 4 β_k’s to reproduce these values. In the second step, we fit the parameters in the momentum-dependent and spin-orbit interactions (ζ and W₀ in [5]) to E/A (energy per nucleon) and R_c (charge radius) of ⁴⁰Ca, ⁴⁸Ca and ²⁰⁸Pb. Agreement to the 6 nuclear data is measured by the χ62 defined by

Among the 7 EoS parameters, we fix ρ₀, E_B, Q_sym to 0.16 fm⁻³, 16.0 MeV and 650 MeV, respectively, and investigate the dependence on the 4 EoS parameters K₀, J, L, and K_sym. In order to avoid too large parameter space, 4 EoS parameters are restricted to the ranges allowed by either experiments or other theories. For K₀, the range is chosen as 220 − 260 MeV, and the interval is divided into 10 MeV, i.e. 220, 230, 240, 250 and 260 MeV. Ranges of J and L are 30 − 34 MeV and 40 − 70 MeV, respectively, and the interval is assumed 1 MeV for both J and L. As for K_sym we consider K_τ

and assume three values −360, −420 and −480 MeV for K_τ . Total number of the combinations of (K₀, J, L, K_τ ) is 5×5×30×3 = 2250, and the resulting χ62 ranges from 7.5 × 10−6 to 3.2 × 10⁻⁴. Among the 2250 results, we plot the results of χ62 less than 10⁻⁴ for each K₀ value in Fig. 1. The results are comparable with the SLy4 and the GSkI models which give χ62 values 2.24 × 10⁻⁴ and 6.9 × 10⁻⁵, respectively.

Figure 1. (Color online) χ62 less than 10⁻⁴ for K₀ = 220, 230, 240, 250 and 260 MeV.

There are several notable features in Fig. 1. First, minimal χ62 is likely to be located in the range K₀ = 230− 250 MeV. This range is consistent with, and narrower than the one from experiment K₀ = 200 − 260 MeV. Second, for larger K₀ values, smaller χ62 is obtained with larger K_τ values (smaller in magnitude) and vice versa. Third, dependence on K_τ is relatively weak in the range K₀ = 230−250 MeV. Interestingly this range agrees with the range of minimum χ62. Since the nuclei in the fitting are not so asymmetric i.e. δ² ≪ 1, and the densities relevant to the 6 data may be ρ₀ and nearby, i.e. j(ρ − ρ₀)/3ρ₀j ≪ 1, it should be natural and reasonable to have weak dependence on K_τ .

From the result of fitting to the nuclear data, one can guess that true K₀ value may be located in K₀ = 230 − 250 MeV. However, it is premature to draw a definite conclusion because density of nuclei is identically close to ρ₀, so the information about the density dependence of the EoS one can obtain from the nuclear data might be limited. In addition, asymmetry δ of ⁴⁸Ca and ²⁰⁸Pb are 0.17 and 0.21, respectively, so they are not large enough to magnify the effect of the symmetry energy. Neutron stars provide environment with more wider density range and the large neutron-proton asymmetry. Not abandoning the possibility of small or large K₀ value, we include K₀ = 220 and 260 MeV in the consideration of the neutron star. For each K₀ value, we select (J, L, K_τ ) that give the two smallest χ62 values. Equation of state parameters thus determined are summarized in Tab. 1. Ranges of J and L in Tab. 1 are 30−34 MeV and 40−70 MeV, respectively. These ranges are similar to the experimental ranges in Ref. [1].

In 1990’s measurement of both mass and radius of a neutron was scarce, so the mass observation played a dominant role in constraining the nuclear EoS at high densities. Maximum mass was known to be about 1.5M_⊙. However masses close to or larger than 2M_⊙ have begun to be reported in the 21st century, and now the general consensus is that the maximum mass of the neutron star may be 2M_⊙ or more. In the core of 2M_⊙ neutron star, density is large enough that the nucleons are overlapping to each other. It is likely that non-nucleonic state will emerge and exist in the core of 2M_⊙ neutron stars, so the EoS at the core is hard to determine accurately because of large uncertainties from various sources.

Table 1 EoS parameters (J, L, K_τ ) giving the two smallest χ62 values for each K₀ value.

K0	(J, L Kτ )	χ2(×10−5)
220	(33, 50, −480)	9.45
	(34, 63, −480)	8.61
230	(33, 66, −420)	3.04
	(33, 52, −480)	3.01
240	(32, 68, −360)	0.75
	(32, 58, −420)	0.89
250	(30, 41, −360)	1.50
	(31, 58, −360)	1.43
260	(30, 47, −360)	5.55
	(31, 63, −360)	6.03

On the other hand, density at the center of a 1.4M_⊙ neutron star does not exceed 3ρ₀ regardless of the stiffness of EoS. At densities smaller than 3ρ₀, uncertainty due to the non-nucleonic degrees of freedom can exist, but their effects are not significant or critical. For example, mass and radius of neutron stars with Λ hyperons were calculated with EoSs over a wide range of stiffness for the nucleon-nucleon (NN), nucleon-Lambda (NΛ) and Lambda-Lambda (ΛΛ) interactions [7]. Maximum mass and the mass-radius curves depend on the NΛ and ΛΛ interactions critically, but the mass-radius relation of the canonical neutron star is seldom affected by the existence of Λ hyperons and the stiffness of their interactions with N and Λ. Therefore 1.4M_⊙ neutron star provides a good laboratory in which we can study the EoS of nucleonic matter at densities higher than ρ₀.

Now more ample and exact data of mass and radius of neutron stars whose mass is close to 1.4M_⊙ are available. We collect the state-of-the-art data of the neutron star in Fig. 2. Horizontal bands around 2M_⊙ exhibit the data of large masses [8–10]. Gray zone represents the range obtained from the low mass X-ray binary (LMXB) data [11], error bars denote the analysis of a soft X-ray source observed in the NICER project [12,13], and the red and the blue regions in the shape of wings correspond to the result of GW170817 (lower small wings) [14] and GW190425 (upper large wings) [15]. In the case of NICER two analyses were performed independently for a single object, so the two results are shown in the figure. A red horizontal line is drawn at mass 1.4M_⊙ for reference. Observation data span broadly on the massradius plane, but there is a narrow band in which all the data overalp. For the mass 1.4M_⊙, the overlapping region restricts the radius in the range 11.8 − 12.5 km.

Figure 2. (Color online) Neutron star mass and radius relations with K₀, J, L and K_τ selected in Tab. 1.

Ten candidates for the EoS parameters (K₀, J, L, K_τ ) also show widely spread behavior. At first compared to the LMXB data at the canonical mass (red line), 4 EoSs are incompatible: 2 soft EoSs corresponding to (K₀, J, L, K_τ ) = (220, 33, 50, −480), (230, 33, 52, −480), and 2 stiff EoSs with (K₀, J, L, K_τ ) = (240, 32, 68, −360), (260, 31, 63, −360). It is notable that K₀ and K_τ are small simultaneously for the 2 soft EoSs, and they are large simultaneously for the 2 stiff EoSs. On the other hand, we have 4 EoSs that fall into the common range of LMXB, NICER and GW for M = 1.4M_⊙: (K₀, J, L, K_τ ) = (230, 33, 66, −420), (240, 32, 58, −420), (250, 31, 58, −360), and (260, 30, 47, −360). There are several interesting features. First, the smallest K₀ and K_τ values, 220 MeV and −480 MeV are ruled out. This provides lower bounds of these parameters. Second, similar to the 2 soft and 2 stiff EoSs incompatible with the LMXB data, smaller K₀ values are grouped with smaller K_τ values, and vice versa. On the other hand, J and L values are small when K₀ and K_τ are large. This means J and L contribute to EoS oppositely to K₀ and K_τ : If J and L soften the EoS, K₀ and K_τ make it stiff. As a result, EoS is adjusted to remain in ranges consistent with data. Most essential aspect of the result is the range of each EoS parameters. If the consistency with neutron star observation is accounted, acceptable ranges could be specified as

all in the unit of MeV. On the other hand, if we take into account the consistency with the nuclear data in addition to the nuclear star observation, K₀ = 260 MeV could be excluded because χ62 is much larger than the other K₀ values. Then we obtain narrower ranges

The ranges are reduced significantly compared to those constrained from experiment [1], and thus they are expected to give EoSs less uncertain from sub to supra saturation densities. Compared to a recent work [16] in which K₀ is fixed to 240 MeV and symmetry energy parameters are fit to 13 nuclear data, J is slightly larger in this work, and L and K_τ ranges are almost identical. With the newly constrained ranges of the EoS parameters, we can test the extreme properties of nuclei such as the neutron skin thickness or the position of the neutron drip line. Consideration about the finite nuclei follows in the next section.

Ranges of the EoS parameters in Eq. (7) are determined using 6 nuclear data of 40,⁴⁸Ca, ²⁰⁸Pb, and the neutron star mass and radius. Since we are aiming at a unified description of both nuclear matter and nuclei, it is mandatory to check the predictive power of our approach for the nuclear properties. Top prior quantity among numerous nuclear properties might be the binding energy.

Table 2 shows the binding energy per nucleon for standard spherical magic nuclei. Numbers in the parentheses denote the difference from experiment in units of %. Energies of 40,⁴⁸Ca and ²⁰⁸Pb are used in fitting the 2 Skyrme force parameters that are coefficients of the grandient and spin-orbit interactions, so it is natural to have good agreement with these data. Energies of ¹⁶O, ⁹⁰Zr and ¹³²Sn are, on the other hand, predictions of the model. Predictions agree well with experiment, giving difference from experiment at the order of 0.5% or less. Another positive aspect of the result is that the predicted energy values are independent of the EoS parameters (K₀, J, L, K_sym). This means that the 4 EoS parameter sets have equal predictive power for the binding energy of magic spherical nuclei.

Table 2 Binding energy per nucleon in MeV for (K₀, J, L, K_τ ) values satsifying the neutron star mass-radius observation. Set A, B, C and D correspond to (230, 33, 66, −420), (240, 32, 58, −420), (250, 31, 58, −360), and (260, 30, 47, −360), respectively. Numbers in parenthesis denote the difference from experiment in %.

	16O	40Ca	48Ca	90Zr	132Sn	208Pb
Exp.	7.976	8.551	8.667	8.710	8.355	7.868
A	7.946 (0.38)	8.564 (−0.16)	8.680 (−0.15)	8.681 (0.33)	8.378 (−0.28)	7.871 (−0.04)
B	7.940 (0.45)	8.555 (−0.04)	8.673 (−0.07)	8.676 (0.39)	8.377 (−0.26)	7.869 (−0.01)
C	7.935 (0.51)	8.546 (0.06)	8.667 (−0.01)	8.678 (0.37)	8.375 (−0.24)	7.868 (0)
D	7.933 (0.54)	8.538 (0.15)	8.653 (0.17)	8.674 (0.42)	8.374 (−0.22)	7.867 (0.02)

One main goal of the work is to have reliable predictions about the objects in extreme conditions in which density is much larger or smaller than the saturation density, or proton-neutron asymmetry approaches to 1. Neutron stars are good laboratory to explore these effects, but nuclei close to the neutron drip line also provide a unique test ground. In Fig. 3, we show the result of two neutron separation energy S_2n of O, Ca, Zr, Sn and Pb isotopes calculated with EoS parameter sets A, B, C, D in Tab. 2. Experimental data are denoted with red squares. In order to account for the pairing correlation consistently in the neutron rich limit, we use the Hartree-Fock-Bogolyubov theory.

Figure 3. (Color online) Two neutron separation energies for the O, Ca, Sn and Pb isotopic chains for the four sets of (K₀, J, L,K_τ ) in Tab. 2.

All the isotopes are assumed spherical, and pairing force is incorporated in the contact interaction form

Parameters in the pairing force are fit to the mean neutron gap of ¹²⁰Sn, 1.392 MeV.

Agreement to the experiment is not as good as the energies of the magic nuclei in Tab. 2. More investigation is necessary for a better reproduction of the data. For example, the disagreement can originate from the way to treat the pairing force. Non-spherical shape of the isotopes could contribute to the discrepancy. Nevertheless, there are a few points that are worthy to be discussed. For the oxigen, results of the 4 EoS parameter sets agree to each other fairly well, but there are sizable differences from experiment at N = 8 and N = 10. Position of the neutron drip line is predicted identically at N = 22. For the calcium, agreement to experiement is reasonable for N less than 30, but it becomes worse for N > 30. For N between 30 and 40, however, the 4 EoS sets predict similar results, and for N > 40, predictions are divided into two branches. Sets with K_τ = −420 MeV give S_2n larger than those of K_τ = −360 MeV sets. Consequently neutron drip lines differ significantly, N = 54 for K_τ = −360 MeV and N = 62 for K_τ = −420 MeV. For the zirconium, agreement to experiment becomes worse for N > 40. Four EoS parameter sets agree to each other well for N less than 82, but the prediction splits into two groups for N > 82. Similar behavior is observed for the tin isotopes. For the lead, the same thing happens at another magic number N = 184. Neutron drip of the Pb isotopes spreads over a huge range of the neutron number, from N = 186 to N = 220.

This work was motivated by the observation that the EoS of nuclear matter which is crucial in understanding the properties of extremely neutron-rich nuclei and the neutron star is not determined precisely enough yet. By reducing the uncertainty in the EoS parameters K₀, J, L, and K_sym, it may be feasible to have more exact knowledge about the state of matter at extreme conditions. KIDS density functional provides a framework adequate for this analysis. Four EoS parameters are determined to best reproduce the binding energy and the charge radius of 40,⁴⁸Ca and ²⁰⁸Pb, and the radius of 1.4M_⊙ neutron stars determined from modern astronomy. As a result we obtain the ranges K₀ ~ 230−250 MeV, J ~ 31−33 MeV, and L ~ 55 − 65 MeV. As for K_τ , we used three values −360, −420 and −480 MeV, and find that −480 MeV makes the EoS too soft to satisfy the neutron star observation. Ranges of J and L are consistent with the ranges in a recent work [16] in which K₀ is fixed to 240 MeV and the nuclear data used in the fitting are chosen differently from this work. EoS parameters thus determined are applied to calculating the properties of nuclei. As for the binding energies of spherical magic nuclei ¹⁶O, ⁹⁰Zr and ¹³²Sn, all the 4 EoS parameter sets consistent with the neutron star data give results agreeing to nuclear data within the errors less than 0.5%. Two-neutron separation energy, on the other hand, shows non-negligible discrepancy with experiment, so it needs more study to understand the origin of the discrepancy and make the theoretical prediction better.

From the analysis of this work, it has been shown that the uncertainties of J and L are reduced to the ranges ±1 MeV and ±5 MeV, respectively.

Then a subsequent question follows: Are these uncertainties narrow enough? or How accurately we should determine the symmetry energy parameters?

The result of two-neutron separation energy shows that the neutron drip line can be sensitive to K_τ .

K_τ is determined by both K_sym and L (Eq. (6)). So called the CSkP (consistent Skyrme parameterization) ranges of K₀ and Q₀ are given as 200 ≤ K0 ≤ 240 MeV and −420 ≤ Q₀ ≤ −360 [1]. Employing these ranges we can approximate Eq. (6) as

If L is determined with errors ±5 MeV, we have about 40−45 MeV uncertainty in K_τ . This uncertainty is similar in magnitude to the interval of K_τ determined in this work −420 − −360 MeV. For a better decision of K_τ , it might be necessary to constrain L in regions narrower than ±5 MeV. Aside from the uncetainty of the EoS parameters, it is important to have more correct description of the exotic nuclei. There are many things to consider, e.g. refined treatment of the pairing force, the effect of deformation, and etc. Some of these subjects are under way.

This work was supported by the Daegu University Research Grant 2018.