1 Introduction

To understand the nature of neutron-rich nucleonic matter has been a major scientific goal in both nuclear physics and astrophysics. The density dependence of nuclear symmetry energy E_sym(ρ) has been a major uncertain part of the equation of state (EOS) of neutron-rich matter especially at high densities, see, e.g., collections in Ref. [1]. Reliable knowledge about the E_sym(ρ) has significant ramifications in answering many interesting questions regarding the structure of rare isotopes and neutron stars, dynamics of heavy-ion collisions and supernova explosions as well as the frequency and strain amplitude of gravitational waves from deformed pulsars and/or cosmic collisions involving neutron stars. During the last two decades, significant efforts have been devoted to exploring the E_sym(ρ) using both terrestrial laboratory experiments [2-12] and astrophysical observations [13-20]. Extensive surveys of the extracted constraints on the E_sym(ρ) around the saturation density ρ₀ indicate that the central values of the E_sym(ρ₀ and its slope L=[3 ρ (∂E_sym/∂ ρ)]ρ0 scatter around 31.6 MeV and 58.9 MeV, respectively [14, 21, 22]. At densities away from ρ₀, however, the E_sym(ρ) remains rather unconstrained especially at supra-saturation densities [12].

Interestingly, recent progresses in another seemingly different field may provide additional information about the density dependence of nuclear symmetry energy. Indeed, theoretical and experimental studies of cold atoms have made impressive progress in recent years, see, e.g., Refs. [23-26] for recent reviews, providing reliable information about the universal EOS (E_UG) of unitary gas (UG) interacting via pairwise s-waves with infinite scattering length but zero effective range. The universal E_UG constrains stringently the EOS of pure neutron matter (PNM) at sub-saturation densities, thus provides possibly additional constraints on the nuclear symmetry energy. In fact, it was recently conjectured that the E_UG provides the lower boundary of the EOS of PNM (E_PNM) [27]. Moreover, using a set of known parameters of symmetric nuclear matter (SNM), and taking zero as an upper bound on the curvature K_sym=[9ρ²(∂²E_sym(ρ)/∂ρ²)]ρ0 of E_sym(ρ) at ρ₀, the authors of Ref. [27] obtained a region of E_sym(ρ₀ - L space that is inconsistent with the unitary gas constraints, excluding many E_sym(ρ) functionals currently actively used in both nuclear physics and astrophysics.

Our original purposes were to examine several issues not clearly addressed in version-1 of Ref. [27]. We notice that some of these issues are now discussed in more detail in its revised version. Nevertheless, it is still useful to provide our results and opinions on some of these issues. The derivation of the excluded region in E_sym(ρ₀ - L space by Ref. [27] relies on several assumptions [28]: the most importance of which is the underlying conjecture that E_PNM(ρ)≥ E_UG(ρ) where E_UG(ρ)=ξ E_F(ρ) with ξ being the Bertsch parameter [23-26] and E_F∝ (ρ/ρ₀)^2/3 is the energy of a non-interacting non-relativistic degenerate Fermi gas of neutrons. This conjecture is not quite the same as stating that PNM has a greater energy than the UG at all densities, since the UG is experimentally accessible only at low densities. This conjecture merely employs an algebraic expression motivated by the result that apparently the low-density neutron gas has a higher energy than the UG, and by the inference that the repulsive nature of the three nucleon (NNN) interaction coupled with the finite values of the range and scattering length of the two-body (NN) s-wave interaction will ensure that the energy of PNM remains above this algebraic expression at higher densities. A further assumption is that possibly attractive higher order NN interactions (p-wave, d-wave, etc.) are not important relative to the repulsive character of three- and higher-body interactions. Additionally, neutrons are assumed to remain non-relativisitc in the density range considered. Given the important ramifications of the findings in Ref. [27], we are motivated to critically examine them adopting the same assumptions. We include the third-order terms in density characterized by the skewness coefficients $J_0=27{\rho }_0^3{\partial }^3{E}_0(\rho ){/}^3{|}_{\rho =\rho }{}_{{}_0}$ and $J_{\text{sym}}=27{\rho }_0^3{\partial }^3{E}_{\text{sym}}(\rho )/\partial {\rho }^3{|}_{\rho ={\rho }_0}$ in expanding the E₀(ρ) and E_sym(ρ), respectively. Particularly, we carefully examine the uncertainties of the curvature K_sym of the symmetry energy and the skewness coefficients J₀ and J_sym, taking into account energy density functionals that are consistent with the PNM EOS derived from microscopic calculations, and examine the effects of those uncertainties on the region of E_sym(ρ₀-L space excluded by the unitary gas constraints.

While Skyrme models consistent with microscopic PNM calculations tend to give K_sym in the range -100 to -200 MeV [32, 37], Relativistic Mean Field (RMF) models consistent with microscopic PNM calculations can give positive values of K_sym [32, 31], reflecting a difference in the form of these two classes of energy density functionals. Indeed, some reputable non-relativistic and relativistic energy density functionals in the literature, see, e.g., reviews in Ref. [29-31, 33], predict positive K_sym values and meet all existing constraints including the EOS of PNM within their known uncertain ranges. For example, the TM2 RMF interaction has K_sym = 50 MeV and passes the PNM test of Ref. [31]. To our best knowledge, while the majority of existing models predict negative values for K_sym, there is no fundamental physics principle excluding positive K_sym values. The current situation clearly calls for more studies on the K_sym especially its experimental constraints. Hopefully, ongoing experiments at several laboratories [34] to extract the isospin dependence of nuclear incompressibility and subsequently the K_sym from giant resonances of neutron-rich nuclei will help settle the issue in the near future. Interestingly, very recently by analyzing comprehensively the relative elliptical flows of neutrons and protons measured by the ASY-EOS and the FOPI-LAND collaborations at GSI using a Quantum Molecular Dynamics (QMD) model [35], the extracted values for the slope and curvature parameters are L=59± 24 MeV and K_sym=88± 372 MeV at the 1σ confidence level. After considering uncertainties of other model parameters including the incompressibility of SNM, neutron-proton effective mass splitting (related to the momentum dependence of the symmetry potential), Pauli blocking and in-medium nucleon-nucleon cross sections, it was concluded that L = 59 ± 24 (exp) ± 16(th) ±10 (sys) MeV and K_sym = 0 ± 370(exp) ± 220(th) ± 150(sys) MeV. To our best knowledge, the latter represents the latest and most accurate constraint on K_sym.

2 The lower boundary of nuclear symmetry energy constrained by the universal EOS of unitary Fermi gas

Within the parabolic approximation for the EOS of isospin asymmetric nuclear matter (ANM) in terms of the energy per nucleon E, i.e., E(ρ,δ)=E₀(ρ)+E_sym(ρ)δ²+O(δ⁴), the symmetry energy E_sym(ρ)= 2^-1[∂²E(ρ,δ)/ ∂ δ²]_{δ =0} can be approximated by

where μ=ρ/ρ₀ is the reduced density and δ=(ρ_n-ρ_p)/ρ is the isospin asymmetry of ANM. Using the conjecture E_PNM(ρ)≥E_UG(ρ) [27] and the EOS of unitary gas

where k_F is the neutron Fermi momentum, the lower boundary of symmetry energy can be obtained from

It is necessary to caution that the above lower boundary of E_sym(μ) is estimated based on the conjecture E_PNM(ρ)≥E_UG(ρ) and the assumption that ξ is a constant in the density range we study. As emphasized in Ref. [27], the conjecture is empirical in nature. While there are strong supports for the conjecture by comparing the EOSs of PNM calculated from various microscopic many-body theories with the E_UG(μ) using a constant ξ≈ 0.37 (see Fig.1 of Ref. [27] and Fig.2 of Ref. [36]) up to about ρ₀, the rigorous condition for unitarity is expected to be reached in PNM only at very low densities. Although one can not prove the validity of the conjecture at high densities, strong physical arguments were made to justify and use it up to about 1.5ρ₀ in Ref. [27]. Thus, the results of our study should be understood with the caveat that they are obtained under the above reasonable but not rigorously proven conjecture and assumptions. Nevertheless, they are useful for comparing with the results of Ref. [27] obtained using the same assumptions.

Figure 1:Full-screen View

(Color online) The skewness parameter J₀ versus the incompressibility K₀ for symmetric matter (a) and the total curvature parameter K_n = K₀ + K_sym versus the total skewness parameter J_n = J₀ + J_sym (b) for all 173 Skyrme and 101 RMF models examined by Dutra et al [30, 31] which pass their pure neutron matter constraints and satisfy 190<K₀<270 MeV.

Figure 2:Full-screen View

(Color online) The lower boundary of symmetry energy as a function of density for different skewness coefficients J₀ = 400, 0, -400, and -800 MeV. The red dashed region represents the variation of E_sym(μ) with K₀ = 230 MeV and J₀ = 0 MeV by adopting ξ = 0.37±0.005. The shadowed region shows the excluded region after considering the uncertainties of ξ, K and J₀.

The EOS of SNM around ρ₀ can be expanded to the third order in density as

in terms of the incompressibility K₀ and skewness J₀. At the saturation point of SNM, we adopt E₀(ρ₀) = -15.9 MeV and ρ₀ = 0.164 fm^-3 [37]. The lower boundary of E_sym(ρ) thus depends on the values of K₀, J₀ and ξ.

The incompressibility K₀ of SNM has been extensively investigated [33, 38], and the most widely used values are K₀ = 240±20 MeV [39, 40] or 230±40 MeV [41]. However, the skewness coefficient J₀ is still poorly known [42-47]. In Fig. 1(a), we show K₀ and J₀ from 274 parameterizations of the Skyrme and RMF models that pass the PNM tests of Dutra et al [30, 31] and satisfy K₀ = 230±40 MeV. The spread in values for J₀ is very large, covering the range ≈ -800<J₀<400 MeV.

As reviewed recently in Refs. [26, 27], currently the best estimate for the Bertsch parameter ξ from lattice Monte Carlo studies is ξ = 0.372(5) consistent with the most accurate experimental value of ξ = 0.376(4). While its values from various other models and experiments have scattered between 0.279 and 0.449(9) within the last decade, it appears that it now has converged to ξ = 0.37±0.005 which we adopt in this work. Shown in Fig. 2 with the red dashed lines are the variation of E_sym(ρ) with ξ = 0.37±0.005, J₀ = 0 and K=230 MeV. Effects of varying the ξ value are very small within the range considered.

Secondly, effects of the skewness parameter are shown by varying the value of J₀ between -800 MeV and 400 MeV, the range covered by the models plotted in Fig. 1(a). Although the range is very large, it translates to a range of uncertainty for E_sym(ρ) that is equivalent to the range of uncertainty in K₀ (80 MeV). This is easy to understand as the expansion of SNM’s EOS converges quickly around the normal density by design (the J₀ contribution of J₀/162 is a factor of 9 less than the K₀/18 term).

Considering the uncertainties of all relevant parameters involved, the most conservative lower boundary of E_sym(ρ) shown as the shadowed region in Fig. 2 is obtained by using ξ=0.37, ρ₀=0.157 fm^-3, E₀(ρ₀)=-15.5 MeV, and K₀=270 MeV; for μ≤ 1, J₀=-800 MeV and for μ>1, J₀=400 MeV. Overall, our observations and results are consistent with the findings in Ref. [27].

3 Constraining the E_sym(ρ₀ versus L boundary

The symmetry energy E_sym(μ) can be expanded around ρ₀ to third order in density as

in terms of its magnitude E_sym(ρ₀, slope L, curvature K_sym and skewness J_sym at ρ₀. Inserting the above equation into Eq. (3), the lower boundary of E_sym(ρ₀ can be expressed as

where K_n=K_sym+K₀ and J_n=J_sym+J₀. Taking the derivative of the above equation with respect to density on both sides, one can readily get an expression for the lower boundary of L,

Then, putting the above expression back to Eq. (6) the latter can be rewritten as

These two equations reveal the correlation between the E_sym(ρ₀ and L along their lower boundaries through the arbitrary density μ. Setting J_n=0, the Eqs. (7) and (8) reduce exactly to the parametric equations of E_sym(ρ₀ and L derived slightly differently in version-1 of Ref. [27]. We note that the quantities that determine the boundary of allowed values of E_sym(ρ₀) and L are the total curvature parameter K_n and total skewness parameter J_n also emphasized in version-2 of Ref. [27].

While having noted that K_sym is experimentally and theoretically poorly known, the E_sym(ρ₀ versus L correlation along their boundaries was obtained in Ref. [27] by setting K_sym=0 based on the prediction of a chiral effective field theory. It was found that the resulting correlation excludes many of the currently actively used models for E_sym(ρ). We reexamine this correlation by varying the ξ, J_n and K_sym within their known uncertain ranges. Again, the value of ξ is now well settled around 0.37± 0.005. Taking K_sym=0, Jn=0 and K₀=230 MeV, the two red dashed lines obtained with ξ = 0.37± 0.005 in both (a) and (b) of Fig. 3 show the resulting lower boundaries of the E_sym(ρ₀ versus L correlation.

Figure 3:Full-screen View

(Color online) The lower boundaries of symmetry energy parameters based on different parameter values. The effects of the uncertainty in J_n (a) and K_sym (b) are demonstrated by the blue lines, holding K_n=230 MeV in (a) and J_n=0 MeV in (b). The red dashed region represents the deviation of E_sym with J₀ = 0 MeV by adopting ξ = 0.37±0.005. The shadowed region shows the excluded region after considering all the uncertainties of ρ₀, E₀(ρ₀), ξ, K₀, J_n, and K_sym(ρ₀).

The skewness coefficients J₀ and J_sym in J_n are both poorly known. We show in Fig 1(b) values of J_n against K_n for the 275 Skyrme and RMF models. J_n varies approximately in the range -500 MeV ≤ J_n ≤ 1000 MeV. To our best knowledge, there is no experimental constraint available on this quantity. K_n varies approximately in the range -150 MeV ≤ K_n ≤ 370 MeV. Most of this comes from the big uncertainties in determining the value of K_sym, which are discussed in detail in Ref. [12]. This is partially because the K_sym depends on not only L but also its derivative (dL/d ρ)_ρ0 by definition. Microscopically, it depends on not only the nucleon isoscalar effective mass $m_0^{*}$ and neutron-proton effective mass slitting $m_{\text{n}}^{*}-m_{\text{p}}^{*}$ but also their momentum and density dependences that are all essentially completely unknown [12]. The latest calculations within many Skyrme Hartree-Fock and/or relativistic mean-field models indicate that -400 ≤ K_sym ≤ 100 MeV [29-31, 33].

The results using K_sym=0, K₀=230 MeV (so K_n = 230 MeV) and ξ=0.37 for different values for J_n are shown in Fig. 3 (a). It is seen that L becomes larger as J₀ decreases, and that the upper boundary of the allowed region will correspond to the lower limit of J_n. At E_sym(ρ₀) = 40 MeV, for example, the increase is only about 20%.

By setting K_sym= -200, -100, 0, and 100 MeV with K₀=230 MeV (corresponding to K_n=30, 130, 230 and 330 MeV) we can illustrate effects of the K_sym on the lower boundary of E_sym(ρ₀ versus L correlation in Fig. 3 (b). It is seen that the K_sym affects the results significantly; as K_sym (and hence K_n) increases, L increases, and the upper boundary of the allowed region will correspond to the largest value of K_sym (and hence K_n). The overall uncertainty in K_n leads to a variation of the upper boundary of the allowed region that is about twice that caused by the uncertainty in J_n. We note that many of the models allowed by K_sym= 0 MeV would be excluded by using K_sym=-200 MeV.

Adopting the following values for the five parameters after taking into account the full range presented by the 275 Skyrme and RMF models, i.e.,

$\begin{array}{l}{\rho }_0=0.157{\text{fm}}^{-3}\mathrm{,}{E}_0({\rho }_0)=-15.5\text{MeV}\mathrm{,}\xi =\mathrm{0.37,}\\ K_{\text{n}}=370\text{MeV}\mathrm{,}{J}_{\text{n}}=-500\text{MeV}\mathrm{,}\end{array}$

a lower boundary excluding only the shadowed region in Fig. 3 is obtained. It is seen that only the TMA and NLρδ, NL3 and LS220 may be excluded, while the STOS, TM1, NLρ, LS220, and KVR, which have been surely excluded previously in Ref. [27] may be allowed. This is mainly a result in extending the upper bound in the uncertainty region of K_sym, since the additional uncertainty in J_n moves the excluded region to the right, (because $J_{\text{n}}\gtrsim 0$). As pointed out in the final analyses of Ref. [27], we agree that strong empirical correlations among J_sym, K_sym and L exist. These correlations can be used to refine the E_sym(ρ₀-L constraint shown in Fig. 3.

It is well known that the detailed density dependence of nuclear symmetry energy E_sym(ρ) contains many interesting and some unknown physics. It is probably not surprising that the estimation of the correlation between its zeroth-order and first-order density expansion coefficients E_sym(ρ₀ and L depends on what we assume about the immediate next high-order term characterized by the curvature K_sym at ρ₀.

4 Concluding remarks

The universal EOS E_UG of the unitary Fermi gas was conjectured in Ref. [27] to provide the lower boundary of the EOS E_PNM of PNM, and thus a constraint on nuclear symmetry energy. Although unproven, the conjecture has strong empirical supports, and its implications are important enough that they should be examined rigorously. We found that the E_UG does provide a useful lower boundary of nuclear symmetry energy. Moreover, this boundary is essentially not affected by the known uncertainty in the skewness coefficient J₀ of SNM. However, it does not tightly constrain the correlation between the magnitude E₀(ρ₀) and slope L unless the curvature K_sym of the E_sym(ρ) at ρ₀ and the skewness parameters J₀ and J_sym are better known. Of these, K_sym is the more important quantity, its uncertainty affecting the lower boundary of E_sym(ρ) by up to twice as much as the uncertainty in Jn. Most of the previously excluded E_sym(ρ) functionals by the universal EOS of unitary Fermi gas assuming K_sym=0 may not be excluded considering the currently known big uncertainties of the K_sym(ρ₀).

For many purposes in both nuclear physics and astrophysics, it is necessary to map out precisely both the E_sym(ρ) and the EOS E₀(ρ) of SNM in a broad density range. Near the saturation density ρ₀, this requires accurate knowledge of the K_sym and J₀ besides the E₀(ρ₀), L and K₀. To this end, it is interesting to mention briefly quantities that are sensitive to the higher-order EOS parameters and current efforts to determine them. For example, the skewness coefficient J₀ characterizes the high-density behavior of E₀(ρ), and it has been found to affect significantly the maximum mass of neutron stars [47]. Moreover, at the crust-core transition point where the incompressibility of neutron star matter at β-equilibrium vanishes, the value of J₀ influences significantly the exact location of the transition point [48]. Thus, astrophysical observations of neutron stars can potentially constrain the J₀ albeit probably not before other EOS parameters are well determined. On the other hand, in terrestrial laboratory experiments, there have been continued efforts to determine the K_sym [33]. One outstanding example is the measurement of the isospin dependence of nuclear incompressibility K(δ)≈K₀+K_τδ²+𝒪(δ⁴) where K_τ=K_sym-6L-J₀L/K₀ using giant resonances of neutron-rich nuclei [49, 50]. While the current estimate of K_τ≈ -550± 100 MeV [33] from analyzing many different kinds of terrestrial experiments is still too rough to constrain tightly the individual values of J₀ and K_sym, new experiments with more neutron-rich beams have the promise of improving significantly the accuracy of the measured K_τ [34]. Thus, we are hopeful that not only the zeroth and first-order parameters K₀, E_sym(ρ₀ and L but also high-order coefficients J₀ and K_sym can be pinned down in the near future by combining new analyses of upcoming astrophysical observations and terrestrial experiments.