Combined Constraints on the Equation of State of Dense Neutron-rich Matter from Terrestrial Nuclear Experiments and Observations of Neutron Stars

First of all, we illustrate the broad variation of E_sym(ρ) by varying independently the coefficients in Equation (3) within their known uncertainty ranges. It is well known that the isospin asymmetry δ(ρ) in neutron stars at a given baryon density ρ is uniquely determined by the E_sym(ρ) in Equation (1) through the charge neutrality and β equilibrium conditions as we shall recall more formally in Section 2.3. Generally speaking, because of the ${E}_{\mathrm{sym}}(\rho )\cdot {\delta }^{2}$ term in the EOS, a higher value of E_sym(ρ) will lead to a smaller δ(ρ) at β equilibrium. As quantitative examples, shown in Figure 1 are the E_sym(ρ) and δ(ρ) as functions of reduced density by varying only one coefficient each time while fixing all others: (a) L = 40, 50, 60, 70, and 80 MeV; (b) K_sym = −400, −300, −200, −100, 0, and 100 MeV; and (c) J_sym = −200, 0, 200, and 400 MeV. As their names indicate, the slope L, curvature K_sym, and skewness J_sym of symmetry energy play different roles and in order become increasingly more important at higher densities. Obviously, variations of them within their currently known uncertainty ranges allow us to sample very different behaviors of the E_sym(ρ) and the corresponding δ(ρ).

Zoom In Zoom Out Reset image size

It is worth noting that some combinations of the parameters lead to a decreasing E_sym(ρ) that may even become negative at high densities. As summarized earlier in Szmaglinski et al. (2006) and reviewed very recently in Li et al. (2018), such a supersoft E_sym(ρ) at high densities was predicted in a number of theoretical calculations using various interactions. At very high densities, when the short-range repulsive tensor force due to the ρ-meson exchange makes the EOS of SNM increase faster with density than that of pure neutron matter where the tensor force is much weaker, the E_sym(ρ) decreases or even becomes negative at high densities (Pandharipande et al. 1972; Wiringa et al. 1988; Li et al. 2008). To the best of our knowledge, such a seemingly unusual high-density behavior of the E_sym(ρ) is not excluded by either any known fundamental physics principle or experiments/observations so far (Tews et al. 2017). Possible consequences of such symmetry energies are discussed in more detail in Szmaglinski et al. (2006), Li et al. (2018), and references therein. In fact, EOSs with a supersoft E_sym(ρ) can still support massive neutron stars if the SNM parts of the EOSs are sufficiently stiff even without the help from the new light mesons proposed and/or possible modified strong-field gravity for massive objects (see, e.g., Krivoruchenko et al. 2009; Wen et al. 2009; Lin et al. 2014; Jiang et al. 2015). Interestingly, while not completely settled yet (Li 2017), there are indeed some circumstantial evidences from intermediate-relativistic energy heavy-ion collisions indicating that the E_sym(ρ) may become supersoft above about 2ρ₀ (Xiao et al. 2009). Currently, devoted efforts are being made by the intermediate-relativistic heavy-ion reaction community to pin down the high-density behavior of nuclear symmetry energy (see, e.g., Li et al. 2014; Russotto et al. 2016; Xu et al. 2016; Trautmann & Wolter 2017; Tsang et al. 2017).

As we shall discuss in detail next, the E_sym(ρ) around the crust–core transition is mostly controlled by the L and K_sym parameters when the E_sym(ρ₀) is fixed at its most probable empirical value of 31.6 MeV. The variation of L from 40 to 80 MeV and of K_sym from −400 to 100 MeV allows us to sample the usual behavior of E_sym(ρ₀) predicted by various nuclear many-body theories in the sub-saturation density region and within their known empirical constraints at ρ₀.

Although the crust of neutron stars contributes only a small fraction of the total mass and radius, it plays an important role in various astrophysical phenomena (Baym et al. 1971a, 1971b; Pethick & Ravenhall 1995; Link et al. 1999; Lattimer & Prakash 2000, 2001, 2007; Steiner et al. 2005; Chamel & Haensel 2008; Sotani et al. 2012; Newton et al. 2013; Pons et al. 2013; Piekarewicz et al. 2014; Horowitz et al. 2015). Critical to many effects of the crust is the transition density ρ_t between the inner crust and the outer core of neutron stars. Previous studies have found consistently that the transition density is very sensitive to the density dependence of nuclear symmetry energy. In particular, the role of the slope parameter L has been extensively studied (see, e.g., Douchin & Haensel 2000; Kubis 2004, 2007; Lattimer & Prakash 2007; Oyamatsu & Iida 2007). Often the studies employ predictions of a particular nuclear many-body theory where the values of L and K_sym are normally correlated. Here we shall study first the individual roles of the L and K_sym in determining the core–crust transition properties and then contours of the transition density and pressure in the L versus K_sym plane. Finally, effects of the L–K_sym correlation based on the systematics from analyzing over 500 nuclear energy density functions are examined.

In the present study, we employ the thermodynamical approach (Kubis 2004, 2007; Lattimer & Prakash 2007) to estimate the crust–core transition point where the uniform matter becomes unstable against being separated into a mixture of single nucleons and their clusters. The method is known to slightly overestimate the transition density compared to the dynamical approach (Xu et al. 2009; Ducoin et al. 2011; Providência et al. 2014), but it is sufficiently good for the purposes of this work. Specifically, the transition density is determined by the vanishing effective incompressibility of neutron star matter at β equilibrium under the charge neutrality condition (Kubis 2004, 2007; Lattimer & Prakash 2007), i.e.,

$$\begin{eqnarray}{K}_{\mu } & = & {\rho }^{2}\displaystyle \frac{{d}^{2}{E}_{0}}{d{\rho }^{2}}+2\rho \displaystyle \frac{{{dE}}_{0}}{d\rho }+{\delta }^{2}\\ & & \times \left[{\rho }^{2}\displaystyle \frac{{d}^{2}{E}_{\mathrm{sym}}}{d{\rho }^{2}}+2\rho \displaystyle \frac{{{dE}}_{\mathrm{sym}}}{d\rho }-2{E}_{\mathrm{sym}}^{-1}{\left(\rho ,\displaystyle \frac{{{dE}}_{\mathrm{sym}}}{d\rho }\right)}^{2}\right]=0.\end{eqnarray} \tag{ 6 }$$

This approach has been used extensively in the literature to locate the transition point using various EOSs (see, e.g., Xu et al. 2009; Ducoin et al. 2011; Providência et al. 2014; Routray et al. 2016). Enclosed in the brackets of the above expression for K_μ are the first-order and second-order derivatives of the symmetry energy, i.e., quantities directly determining the L and K_sym. It is thus necessary and interesting to first explore separate roles of the latter on the transition density. Shown in panels (a) and (b) of Figure 2 are the transition density ρ_t as functions of L with different values of K_sym and K_sym with different values of L, respectively. It is clearly seen that the ρ_t changes more dramatically with the variation of K_sym than L in their respective uncertainty ranges. This is mainly because the last two terms in the expression for K_μ largely cancel out, leaving the curvature of E_sym(ρ) dominant. In addition, the value of L is already relatively well constrained in a smaller range than the K_sym, making the variation of ρ_t with L look weaker.

Zoom In Zoom Out Reset image size

Next, we examine the transition density ρ_t and pressure P_t by varying both the K_sym and L within their uncertainty ranges continuously. Shown in the two panels of Figure 3 are contours of constant transition densities ρ_t and pressures P_t in the L–K_sym plane, respectively. For transition densities larger than ρ_t = 0.07 fm⁻³, the required K_sym increases monotonically with L, while different behaviors are observed for lower values of ρ_t. The lowest transition density about ρ_t = 0.0549 fm⁻³ appears around the boundary corner at L = 77 MeV and K_sym = 100 MeV. Different from the contours of constant ρ_t, for a fixed P_t, the required K_sym always increases linearly with L before reaching the P_t = 0 boundary along the line K_sym = 3.64 L −163.96 (MeV). The latter is used as a limit in exploring properties of neutron stars in the EOS parameter space.

Zoom In Zoom Out Reset image size

As indicated earlier, in using Equations (2) and (3) we assume that the coefficients are independent and intend to constrain them directly from observations using as little as possible predictions of any particular many-body theory. Thus, in determining the crust–core transition point for constructing the EOS of neutron stars, we freely vary the K_sym and L within their uncertainty ranges specified earlier. Nevertheless, theoretically predicted values of the K_sym and L are often correlated when model ingredients and/or interactions are varied. Thus, for a consistency check we also study how the predicted correlation between K_sym and L may limit the transition point. As discussed by many people in the literature (see, e.g., Farine et al. 1978; Chen et al. 2009; Danielewicz & Lee 2009; Vidaña et al. 2009; Ducoin et al. 2011; Pearson et al. 2014; Providência et al. 2014; Mondal et al. 2017; Tews et al. 2017), based on the systematics of many predictions using various many-body theories and interactions, an approximately universal and linear correlation exists between the K_sym and L. For example, using a total of over 500 energy density functions, including 263 relativistic mean field (RMF) models and Hartree–Fok calculations using 240 Skyrme (Dutra et al. 2012, 2014) and some realistic and Gogny interactions, Mondal et al. (2017) found the following K_sym–L correlation:

$$\begin{eqnarray}&&{K}_{\mathrm{sym}}=(-4.97\pm 0.07)(3{E}_{\mathrm{sym}}({\rho }_{0})-L)+66.8\\ &&\pm 2.14\,\,\mathrm{MeV}.\end{eqnarray} \tag{ 7 }$$

Using essentially the same sets of energy density functionals but requiring 0.149 fm⁻³ < ρ₀ < 0.17 fm⁻³, −17 MeV < E₀(ρ₀) < −15 MeV, 25 MeV < E_sym(ρ₀) < 36 MeV, and 180 MeV < K₀ < 275 MeV, Tews et al. (2017) rejected some of the energy functionals. They found the following K_sym–L correlation using the remaining 188 Skyrme and 73 RMF interactions:

$$\begin{eqnarray}&&{K}_{\mathrm{sym}}=3.501L-305.67\pm 24.26\,(\pm 56.59)\,\,\mathrm{MeV},\end{eqnarray} \tag{ 8 }$$

where the ±24.26 and ±56.59 are error bars including 68.3% and 95.4% of the accepted EOSs, respectively. In Figure 3, the above two K_sym–L correlation functions are shown with the white and red dashed lines, separately. The 68.3% uncertainty range (±24.26 MeV) in Equation (8) is indicated by the red dotted lines, while that of Equation (7) is too small to be shown here. While the parameterizations in Equations (7) and (8) are consistent, their mean values have slightly different slopes because of the different selection criteria used. Nevertheless, if the uncertainty range of Equation (8) is enlarged from 68.3% to 95.4% of the accepted EOSs, the parameterization in Equation (7) can then be fully covered by that in Equation (8). Using the above two correlations, the transition density and pressure are then restricted to be about ρ_t = 0.08 fm⁻³ and P_t = 0.40 MeV fm⁻³, respectively. These values are consistent with the crust–core transition properties often used in the literature. To this end, especially since some of the apparent correlations among EOS parameters from model calculations may be spurious (Margueron et al. 2017a), it is worth noting that while the K_sym–L correlation from the systematics of over 500 energy density functions is very useful for the consistency check, it is still necessary and important to determine the individual values of K_sym and L from experiments/observations. In our following calculations, we thus consistently use the crust–core transition density and pressure by varying independently the K_sym and L values within their respective uncertain ranges without using any of the above two correlation functions.

For completeness and the ease of our discussions in the following, we first recall here the formalism for calculating the EOS in the cores of neutron stars. The total energy density (ρ, δ) of charge-neutral npeμ matter at β equilibrium can be written as

$$\begin{eqnarray}&&\epsilon (\rho ,\delta )={\epsilon }_{b}(\rho ,\delta )+{\epsilon }_{l}(\rho ,\delta ),\end{eqnarray} \tag{ 9 }$$

where _b(ρ, δ) and _l(ρ, δ) are the energy density of baryons and leptons, respectively. The _b(ρ, δ) can be calculated from

$$\begin{eqnarray}&&{\epsilon }_{b}(\rho ,\delta )=\rho E(\rho ,\delta )+\rho {M}_{N},\end{eqnarray} \tag{ 10 }$$

where the specific energy E(ρ, δ) of baryons is given in Equation (1) and M_N is the average rest mass of nucleons. The _l(ρ, δ) from the noninteracting Fermi gas model can be expressed as (${\hslash }=c=1$) (Oppenheimer & Volkoff 1939)

$$\begin{eqnarray}&&{\epsilon }_{l}(\rho ,\delta )=\eta \phi (t)\end{eqnarray} \tag{ 11 }$$

with

$$\begin{eqnarray}&&\eta =\displaystyle \frac{{m}_{l}^{4}}{8{\pi }^{2}},\,\,\phi (t)=t\sqrt{1+{t}^{2}}(1+2{t}^{2})-{ln}(t+\sqrt{1+{t}^{2}}),\end{eqnarray} \tag{ 12 }$$

and

$$\begin{eqnarray}&&t=\displaystyle \frac{{(3{\pi }^{2}{\rho }_{l})}^{1/3}}{{m}_{l}}.\end{eqnarray} \tag{ 13 }$$

The chemical potential of particle i can be calculated from

$$\begin{eqnarray}&&{\mu }_{i}=\displaystyle \frac{\partial \epsilon (\rho ,\delta )}{\partial {\rho }_{i}}.\end{eqnarray} \tag{ 14 }$$

The isospin asymmetry δ(ρ) and relative particle fractions at different densities in neutron stars are obtained through the β equilibrium condition ${\mu }_{n}-{\mu }_{p}={\mu }_{e}={\mu }_{\mu }\approx 4\delta {E}_{\mathrm{sym}}(\rho )$ and the charge neutrality condition ρ_p = ρ_e + ρ_μ. The pressure of the system can be calculated numerically from

$$\begin{eqnarray}&&P(\rho ,\delta )={\rho }^{2}\displaystyle \frac{d\epsilon (\rho ,\delta )/\rho }{d\rho }.\end{eqnarray} \tag{ 15 }$$

The above expressions allow us to calculate the core EOS, which is connected smoothly at the core–crust transition point to the NV EOS (Negele & Vautherin 1973) for the inner crust followed by the BPS EOS (Baym et al. 1971b) for the outer crust.

As mentioned earlier, the EOS parameters K₀, E_sym(ρ₀), and L near ρ₀ are relatively well determined. To investigate how/what high-density EOS parameters are constrained by the three astrophysical observations considered in this work, we construct the EOS of neutron star matter by varying the poorly known J₀, K_sym, and J_sym characterizing the EOS of dense neutron-rich nucleonic matter. In principle, all coefficients used in Equations (2) and (3) should be varied simultaneously within a multivariant Bayesian inference (Steiner et al. 2010; Raithel et al. 2016; Margueron et al. 2017b). Such a study is in progress. In the present work, we shall perform traditional analyses first in the 3D parameter space spanned by J₀, K_sym, and J_sym while fixing all other parameters at their currently known most probable values. Equivalent to re-parameterizing the EOS of SNM and E_sym(ρ) with fewer parameters as often done in the literature, or expanding Equations (2) and (3) only up to ${[(\rho -{\rho }_{0})/3{\rho }_{0}]}^{2}$, we shall also explore the EOS in the 2D parameter space of L and K_sym by setting J₀ = J_sym = 0 while keeping other parameters at their known most probable values. The two cases studied here are similar in spirit to using different numbers of piecewise polytropes or parameters to model the EOS of dense neutron-rich matter. Naturally, the values of the parameters involved may be different in the two cases, but they should all asymptotically approach the same existing constraints on them near ρ₀.

Certainly, there have been continuous efforts in both astrophysics and nuclear physics to constrain the EOS parameters in both cases. For example, from the pressure of SNM constrained by nucleon collective flow data in relativistic heavy-ion collisions (Danielewicz et al. 2002), a constraint of −1280 MeV ≤ J₀ ≤ −10 MeV was obtained by Cai & Chen (2017). Combining it with the mass of neutron star PSR J0348+0432, they further narrowed it down to −494 MeV ≤ J₀ ≤ −10 MeV. By analyzing X-ray bursts, Steiner et al. (2010) extracted a value of −690 MeV ≤ J₀ ≤ −208 MeV or −790 MeV ≤ J₀ ≤ −330 MeV assuming a photospheric radius of r_ph ≫ R or r_ph = R, respectively. All of these constraints on J₀ overlap but have different uncertainty ranges. Similarly, the K_sym and J_sym have not been well constrained either by any experiments/observations so far. Nevertheless, the systematics of over 500 RMF and SHF energy density functionals indicates the following range: −800 MeV ≤ J₀ ≤ 400 MeV (Dutra et al. 2012, 2014), −400 MeV ≤ K_sym ≤ 100 MeV, and −200 MeV ≤ J_sym ≤ 800 MeV, respectively (see, e.g., Chen et al. 2009; Dutra et al. 2012, 2014; Colò et al. 2014; Zhang et al. 2017). Thus, we adopt these ranges for the parameters J₀, K_sym, and J_sym to be consistent with both existing experimental and theoretical findings.

Within the above uncertainty ranges of the EOS parameters, the pressure in neutron stars can be varied within the shaded band shown in the left panel of Figure 4. Its upper and lower limits are obtained by using the parameter set of L = 80 MeV, K_sym = 100 MeV, J_sym = 800 MeV, and J₀ = 400 MeV and the set of L = 40 MeV, K_sym = −400 MeV, J_sym = 134 MeV, and J₀ = 400 MeV, respectively. The individual roles of these parameters are examined by varying them independently in the right panel. As one expects, the variation of J₀ is most effective in modifying the pressure at supra-saturation densities.

Zoom In Zoom Out Reset image size

To be clear, we first emphasize again that in this case the following parameters are fixed at their currently known most probable values based on previous systematic surveys as discussed earlier: K₀ = 240 MeV, E_sym(ρ₀) = 31.7 MeV, and L = 58.9 MeV. We explore constraints on the EOS in the 3D K_sym–J_sym–J₀ parameter space with the crust–core transition density consistently determined and the condition that P_t ≥ 0. Technically, in exploring the 3D parameter space in three loops, we change the J₀ to reproduce specific observational data at given values of K_sym and J_sym that are varied independently. Shown in Figure 5 are the constant surfaces of neutron stars' maximum mass of M_max = 2.01 M_⊙ (green), radius of R_1.4 = 12.83 km (magenta) and 10.62 km (yellow), and the upper limit of the dimensionless tidal deformability Λ_1.4 = 800 (orange) of canonical neutron stars, respectively. For clarity, a causality surface limiting M ≤ 2.4 M_⊙ is not shown.

Zoom In Zoom Out Reset image size

For ease of the following discussions, it is worth recalling first that, as shown in Equation (1), the EOS of SNM E₀(ρ), the symmetry energy E_sym(ρ), and the isospin asymmetry δ(ρ) at β equilibrium are the three quantities together determining the total pressure in neutron stars. More specifically, the total pressure is proportional to ${{dE}}_{0}(\rho )/d\rho +{\delta }^{2}\cdot {{dE}}_{\mathrm{sym}}(\rho )/d\rho$. While the dE₀(ρ)/dρ term is controlled by the J₀ parameter with a fixed K₀, the dE_sym(ρ)/dρ and δ(ρ) are determined by the K_sym and J_sym parameters when L is fixed. Moreover, the symmetry energy contribution to the total pressure is weighted by δ²(ρ). When the E_sym(ρ) is softer with smaller or negative K_sym and J_sym values, the system is more neutron-rich as shown in Figure 1. In particular, for extremely small K_sym and J_sym values (e.g., K_sym = −400 MeV and J_sym = −200 MeV), the E_sym(ρ) becomes negative and the δ(ρ) reaches its maximum of 1 at high densities. Then, the necessary contribution to the pressure from the E₀(ρ) term will require a large J₀ value to support massive neutron stars. At the other extreme, however, when both the K_sym and J_sym are strongly positive (e.g., K_sym = 100 MeV and J_sym = 800 MeV at the bottom left corner), the symmetry energy E_sym(ρ) is superstiff and the δ(ρ) is very small as shown in Figure 1. The required J₀ to support massive neutron stars is then very small.

It is interesting to note several major features in this rather information-rich 3D plot summarizing very extensive calculations. Let us first focus on the constant surface of M_max = 2.01 M_⊙. From the top right to the bottom left corner, the required J₀ first decreases quickly and then remains almost constant with the increasing values of K_sym and J_sym from negative to positive. This feature is completely understandable based on the discussions in the previous paragraph. Namely, near the upper right corner, the symmetry energy is supersoft and the resulting δ(ρ) is close to 1. The weight δ²(ρ) of the symmetry energy contribution to the pressure is significant, while the dE_sym(ρ)/dρ value is small and may even be largely negative. To support neutron stars with the maximum mass of M = 2.01 M_⊙, the value of J₀ has to be highly positive to compensate the small contribution or overcome the possibly negative contribution from the symmetry energy. However, near the bottom left corner where the symmetry energy is superstiff, the resulting δ² at β equilibrium becomes so small that the symmetry energy contribution to the pressure is strongly suppressed. The necessary value of J₀ is therefore small, and the constant surface of M_max = 2.01 M_⊙ becomes rather flat at small J₀ values. More quantitatively, the required minimum value of J₀ is about −243 MeV at K_sym = −400 MeV and J_sym = 800 MeV.

For a comparison and to see more clearly the role of J₀ in determining the masses of neutron stars, a constant surface of M_max = 2.74 M_⊙ (sky blue) corresponding to the composite mass of GW170817 is also shown. It is seen that the overall shapes of the constant surfaces of M_max = 2.01 and 2.74 M_⊙ are rather similar. To support such a hypothetical neutron star would require a J₀ value above 400 MeV beyond its current uncertain range and the causality limit. While the fate of GW170817 is not completely determined yet, several analyses using observations of GW170817 combined with theories/simulations and the causality limit under some caveats have placed the upper bounds on neutron star masses in the range of 2.16–2.28 M_⊙ (Margalit & Metzger 2017; Shibata et al. 2017; Rezzolla et al. 2018; Ruiz et al. 2018; see also Radice et al. 2018, for a very recent review). For instance, Margalit & Metzger (2017) used electromagnetic constraints on the remnant imposed by the kilonova observations after the merger and the gravitational wave information predicted a maximum mass of M_max ≤ 2.17 M_⊙ with 90% confidence. To support neutron stars with such masses, the J₀ just needs to be slightly positive well within its current uncertain range. Thus, once the maximum mass is pinned down, it will put a stringent upper limit on the J₀ parameter.

In addition, it is also known that quadrupole deformations of neutron stars depend on the symmetry energy (Krastev et al. 2008a, 2008b; Fattoyev et al. 2013, 2014, 2017; Krastev & Li 2018). Interestingly, Abbott et al. (2017) inferred that the dimensionless tidal deformability of GW170817 has an upper limit of Λ_1.4 ≤ 800 at the 90% confidence level for the low-spin prior. However, it is seen in Figure 5 that the constant surface of Λ_1.4 = 800 (orange) locates far outside the constant surface of R_1.4 = 12.83 km. Thus, limits on the high-density EOS parameters from the Λ_1.4 ≤ 800 constraint alone are currently much looser than the radius constraint extracted from analyzing the X-ray data. Nevertheless, the expected detection of gravitational waves from a large number of neutron star mergers has the potential to improve the situation.

Now, let us move to the EOS constraints from the radii of neutron stars. The radius of canonical neutron stars is known to depend strongly (weakly) on the nuclear symmetry energy (EOS of SNM; Li & Steiner 2006). It is thus not surprising that the two constant surfaces of radius at R_1.4 = 10.62 km and R_1.4 = 12.83 km for canonical neutron stars are essentially vertical in Figure 5, indicating a weak dependence on the J₀ as one expects. Indeed, they have significant dependences on both the K_sym and J_sym as indicated by the separation between the two constant-radius surfaces. More quantitatively, the required values of J₀ in the constant surfaces of R_1.4 = 10.62 km and R_1.4 = 12.83 km decrease continuously with increasing K_sym and J_sym. The two surfaces can be approximately described by J₀ = −585.64 (MeV) − 2.86 K_sym − 1.00 J_sym and J₀ = 182.54 (MeV) − 3.19 K_sym − 0.60 J_sym, respectively. With a fixed value of L, the nuclear pressure becomes stronger with increasing values of K_sym and J_sym. Thus, the constant surface of R_1.4 = 12.83 km is on the left (having stiffer symmetry energies) of the R_1.4 = 10.62 km surface.

As indicated by the arrows in Figure 5, the space enclosed by the three constant surfaces of M_max ≥ 2.01 M_⊙ and 10.62 km ≤ R_1.4 ≤ 12.83 km is the EOS parameter space allowed by the astrophysical observations of the maximum mass and radii of neutron stars. The space is constrained from the bottom by the maximum mass, and its intersections with the two limits on the radii set the boundary lines restricting the K_sym and J_sym at small J₀ values around −200 MeV. At higher J₀ values, the EOS parameter space is mainly bounded by the two radius constraints. To show the accepted EOS parameters more clearly, the allowed regions in the K_sym versus J_sym planes with J₀ = 400, 200, 0, and −200 MeV, respectively, are shown in Figure 6. The boundaries of these allowed regions in the K_sym–J_sym planes are shown in Figure 7. The regions inside the lines are allowed for a given J₀. The shadowed regions are excluded values of K_sym and J_sym for all J₀ values. More specifically, the boundary of the right shadowed region can be divided into two parts based on the slopes. The upper one with a negative slope is obtained by requiring the crust–core transition pressure to always stay positive, i.e., P_t ≥ 0. It can be described approximately by the expression J_sym = −11.00 K_sym + 457.71 MeV. The lower one with a positive slope is obtained by the intersection line between the surfaces of M_max = 2.01 M_⊙ and R_1.4 = 12.83 km. It can be fitted by the expression J_sym = 7.68 K_sym − 504.90 MeV. This boundary sets an upper limit for K_sym at about 68 MeV. The left shadowed region is excluded by the intersection line of R_1.4 = 10.62 km and J₀ = 400 MeV in Figure 5. It can be fitted by the expression J_sym = −3.07 K_sym − 1054.89 MeV.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Overall, within the framework of our analyses, using the maximum mass of neutron stars and the upper and lower limits of the radii of canonical neutron stars, the three EOS parameters are only limited in a space shown in Figure 5 not completely closed, with its boundaries partially given in Figure 7. Obviously, data of more independent observables from either or both astrophysics and nuclear physics are needed to determine the individual values of K_sym, J_sym, and J₀.

Within the currently known uncertainty ranges of J₀ and J_sym, one can parameterize the EOS with fewer parameters by setting J₀ = J_sym = 0 in Equations (2) and (3) as often done in the literature. Then, the two most poorly known parameters are K_sym and L. The latter is known to be around L ≈ 58.7 ± 28.1 MeV as mentioned earlier. In this 2D model framework, here we explore how/what the same three astrophysical constraints may teach us about the EOS of dense neutron-rich matter.

Shown in Figure 8 are contours of constant maximum masses of neutron stars in the L–K_sym plane. The red, magenta, blue, and cerulean shadows represent regions where M_max ≤ 2.01 M_⊙, R_1.4 ≥ 12.83, R_1.4 ≤ 10.62 km, and Λ_1.4 ≥ 800, respectively. In the white excluded region, the crust–core transition pressure P_t ≤ 0. It is seen that the maximum mass of neutron stars increases with increasing K_sym as one expects. However, it is rather insensitive to the variation of L in the region considered. Again and obviously, the tidal polarizability Λ_1.4 ≤ 800 is much less restrictive than the radius constraint of R_1.4 ≤ 12.83 km. It is also seen that the R_1.4 ≥ 10.62 km constraint is covered by the constraint M_max ≥ 2.01 M_⊙. Consequently, the area bounded by the curves of M_max ≥ 2.01 M_⊙ and R_1.4 ≤ 12.83 km is the region allowed by the astrophysical observations considered. In this region, the maximum value of K_sym is about 52 MeV, consistent with that found in the 3D analyses in the previous subsection. Also in this region, the upper limit of the maximum mass is about 2.37 M_⊙, reached at L ≈ 60 MeV. Probably incidentally, the latter is also currently the most probable value of L based on the surveys of 53 analyses of existing data (Li & Han 2013; Oertel et al. 2017). Interestingly, the observed upper limit of the maximum mass is in good agreement with the findings by Fryer et al. (2015), Lawrence et al. (2015), and Alsing et al. (2017). They found that the upper limit of the maximum mass of neutron stars is between 2.0 and 2.5 M_⊙. However, it is necessary to caution here that we have taken J₀ = J_sym = 0 in the 2D study. As we have discussed in the previous subsection, the value of J₀ affects significantly the maximum mass of neutron stars. Thus, without more precise knowledge about the J₀ parameter, the absolute maximum mass of neutron stars cannot be pinned down. Based on our 2D model analyses here, neutron stars more massive than 2.37 M_⊙ would require a positive value of J₀.

Zoom In Zoom Out Reset image size