Ruling out the supersoft high-density symmetry energy from the discovery of PSR J0740+6620 with mass $2.14^{+0.10}_{-0.09}M_\odot$

TL;DR

This paper uses recent neutron star mass measurements to argue that the high-density symmetry energy cannot be supersoft, challenging previous assumptions and constraining neutron star properties such as tidal deformability.

Contribution

It demonstrates that the recent mass of PSR J0740+6620 rules out supersoft symmetry energy at high densities, providing new constraints on neutron star equations of state.

Findings

Supersoft symmetry energy becomes negative at high densities in neutron stars.

The existence of a 2.14 solar mass neutron star constrains the tidal deformability to be at least 348.

Supersoft high-density symmetry energy cannot explain the observed massive neutron stars.

Abstract

Using the very recently reported mass $2.14^{+0.10}_{-0.09}M_\odot$ of PSR J0740+6620 together with the data of finite nuclei and the constraints on the equation of state of symmetric nuclear matter at suprasaturation densities from flow data in heavy-ion collisions, we show that the symmetry energy $E_{\rm sym}(n)$ cannot be supersoft so that it becomes negative at suprasaturation densities in neutron stars (NSs) and thus may make the NS have a pure neutron matter core. This is in contrast to the fact that using the mass $2.01 \pm 0.04 M_\odot$ of PSR J0348+0432 as the NS maximum mass cannot rule out the supersoft high-density $E_{\rm sym}(n)$. Furthermore, we find the stiffer high-density $E_{\rm sym}(n)$ based on the existence of $2.14M_\odot$ NSs leads to a strong constraint of $\Lambda_{1.4} \ge 348^{+88}_{-51}$ for the dimensionless tidal deformability of the canonical $1.4M_\odot$ NS.

Tables

Table 1. Table 1: Experimental data on the binding energies E B subscript 𝐸 B E_{\rm{B}} (12 spherical even-even nuclei)  (Wang et al., 2017 ) , the charge r.m.s. radii r c subscript 𝑟 c r_{\rm{c}} (9 nuclei)  (Angeli & Marinova, 2013 ; Fricke et al., 1995 ; Le Blanc et al., 2005 ) , the isoscalar giant monopole resonance (GMR) energies E GMR subscript 𝐸 GMR E_{\rm{GMR}} and its experimental error (4 nuclei)  (Youngblood et al., 1999 ) , and 5 spin-orbit energy level splittings ϵ ls A superscript subscript italic-ϵ ls 𝐴 \epsilon_{\rm{ls}}^{A}   (Vautherin & Brink, 1972 ) . Here ν ​ ( π ) 𝜈 𝜋 \nu(\pi) denotes neutron(proton).

AX	$E_B$(MeV)	$r_c$(fm)	$E_{\mathrm{GMR}}$(MeV)	${ϵ}_{\mathrm{ls}}^A$(MeV)
16O	$-127.619$	2.6991		6.10(1p$\nu$)
				6.30(1p$\pi$)
40Ca	$-342.052$	3.4776
48Ca	$-416.001$	3.4771
56Ni	$-483.995$	3.7760
68Ni	$-590.408$
88Sr	$-768.468$	4.2240
90Zr	$-783.898$	4.2694	17.81$\pm$0.35
100Sn	$-825.300$
116Sn	$-988.681$	4.6250	15.90$\pm$0.07
132Sn	$-1102.84$
144Sm	$-1195.73$	4.9524	15.25$\pm$0.11
208Pb	$-1636.43$	5.5012	14.18$\pm$0.11	1.32(2d$\pi$)
				0.89(3p$\nu$)
				1.77(2f$\nu$)

Equations

$$E(n,\delta )=E_0(n)+E_{\mathrm{sym}}(n){\delta }^2+\cdots ,$$

$$E_{\mathrm{sym}}(n)={\frac{1}{2!}\frac{{\partial }^{2}E(n,\delta )}{\partial {\delta }^2}|}_{\delta =0}.$$

$$E_0(n)=E_0(n_0)+\frac{K_0}{2!}{\chi }^2+\frac{J_0}{3!}{\chi }^3+\cdots ,$$

$$E_{\mathrm{sym}}(n)$$

$$v({\mathbf{r}}_i,{\mathbf{r}}_j)$$

$${\Lambda }_M=\frac{2}{3}{k}_2(R/M{)}^5.$$

$$k_2$$

$$\frac{dy(r)}{dr}=-\frac{y(r{)}^2+y(r)F(r)+r^{2}Q(r)}{r},$$

$$F(r)$$

$$Q(r)$$

$$\frac{dP(r)}{dr}$$

$$\frac{dM(r)}{dr}$$

$${\chi }^2(\mathbf{p})=\sum _{i=1}^N{\left(\frac{{\mathcal{O}}_i^{\mathrm{th}}(\mathbf{p})-{\mathcal{O}}_i^{\exp }}{△{\mathcal{O}}_i}\right)}^2,$$

Full text

dense matter — equation of state — stars: neutron — tidal deformability

1 Introduction

The density dependence of nuclear symmetry energy $E_{\rm sym}(n)$, which characterizes the isospin dependence of the equation of state (EOS) of nuclear matter, is fundamentally important due to its multifaceted roles in nuclear physics and astrophysics (Danielewicz et al., 2002; Lattimer and Prakash, 2004; Steiner et al., 2005; Baran et al., 2005; Li et al., 2008; Baldo & Burgio, 2016; Özel & Freire, 2016; Lattimer and Prakash, 2016; Watts et al., 2016; Oertel et al., 2017; Wolter, 2018; Blaschke & Chamel, 2018; Li et al., 2019). Theoretically, it is still a big challenge to calculate the $E_{\rm sym}(n)$ directly from the first-principle non-perturbative QCD (Brambilla et al., 2014), and currently information on the $E_{\rm sym}(n)$ is mainly obtained in the effective models. So far essentially all available nuclear effective models have been used to calculate the $E_{\rm sym}(n)$, and the results can be roughly classified equally into two groups (see, e.g., Refs. (Stone et al., 2003; Chen, 2017)), i.e., a group where the $E_{\rm sym}(n)$ increases with the density $n$, and the other group where the $E_{\rm sym}(n)$ first increases with $n$ and then decreases above a certain suprasaturation density and even becomes negative at high densities. The $E_{\rm sym}(n)$ in the latter group is generally regarded as soft, and here we regard the $E_{\rm sym}(n)$ as supersoft if it becomes negative at the suprasaturation densities inside neutron stars (NSs). In this sense, the supersoft $E_{\rm sym}(n)$ may make the NS have a pure neutron matter (PNM) core, which will have important implications on the chemical composition and cooling mechanisms of protoneutron stars (Lattimer et al., 1991; Sumiyoshi & Toki, 1994; Prakash et al., 1997), the critical densities for the appearance of hyperons (Providência et al., 2019) and antikaon condensates (Lee, 1996; Kubis & Kutschera, 1999) in NSs, the NS mass-radius relations (Prakash et al., 1988; Engvik et al., 1994), and the possibility of a mixed quark-hadron phase (Kutschera & Niemiec, 2000; Wu & Shen, 2019) in NSs.

Unfortunately, the high-density behavior of the $E_{\rm sym}(n)$ is still very elusive, although the $E_{\rm sym}(n)$ at subsaturation densities has been relatively well determined from analyzing the data of finite nuclei (see, e.g., Refs. (Zhang & Chen, 2013; Brown, 2013; Danielewicz & Lee, 2014; Zhang & Chen, 2015). In terrestrial laboratories, the high-density nuclear matter can be produced only in high-energy heavy-ion collisions, and presently the resulting high-density $E_{\rm sym}(n)$ from analyzing the data in heavy-ion collisions can be either supersoft or stiff, strongly depending on the models and data (Xiao et al., 2009; Feng & Jin, 2010; Russotto et al., 2011; Xie et al., 2013; Cozma et al., 2013; Hong & Danielewicz, 2014; Russotto et al., 2016; Zhang & Ko, 2017). In nature, the NSs provide an ideal astrophysical site to explore the high-density $E_{\rm sym}(n)$. In particular, the dimensionless tidal deformability $\Lambda_{M}$ for a NS with mass $M$, which is specifically sensitive to the NS radius and thus the high-density $E_{\rm sym}(n)$, can be extracted from the gravitational wave (GW) signal of the binary neutron star (BNS) merger (Hinderer, 2008; Flanagan & Hinderer, 2008; Hinderer et al., 2010; Vines et al., 2011; Damour et al., 2012). Actually, the limit of $\Lambda_{1.4}\leq 580$ from the recent GW signal GW170817 (Abbott et al., 2017a, b) already can exclude too stiff high-density $E_{\rm sym}(n)$ (Zhou et al., 2019). In addition, the existence of large mass NSs may set a lower limit for the high-density $E_{\rm sym}(n)$, but we note that the observed mass $2.01\pm 0.04M_{\odot}$ of PSR J0348+0432 (Antoniadis et al., 2013) is still consistent with the supersoft high-density $E_{\rm sym}(n)$ (Zhou et al., 2019).

Very recently, a millisecond pulsar J0740+6620 with mass $2.14^{+0.10}_{-0.09}M_{\odot}$ ($68.3\%$ credibility interval) was reported (Cromartie et al., 2019) by combining the relativistic Shapiro delay data taken over 12.5 years at the North American Nanohertz Observatory for Gravitational Waves with recent orbital-phase-specific observations using the Green Bank Telescope. This pulsar may hence replace the previously reported heaviest PSR J0348+0432 with mass $2.01\pm 0.04M_{\odot}$ (Antoniadis et al., 2013) and set a new record for the maximum mass of NSs. It is thus interesting to examine whether this new heaviest NS can give new insight on the high-density $E_{\rm sym}(n)$.

In this work, using the data of finite nuclei together with the constraints on the EOS of symmetric nuclear matter (SNM) at suprasaturation densities from heavy-ion collisions, we show the existence of NSs with mass $2.14M_{\odot}$ can rule out the supersoft $E_{\rm sym}(n)$, although the largest NS mass $2.01M_{\odot}$ cannot. We further find the stiffer lower limit of the high-density $E_{\rm sym}(n)$ from the existence of NSs with mass $2.14M_{\odot}$ leads to a quite large lower bound value for $\Lambda_{1.4}$, i.e., $\Lambda_{1.4}\geq 348^{+88}_{-51}$.

2 Model and method

2.1 Nuclear matter EOS

For an isospin asymmetric nuclear matter with neutron (proton) number density $n_{\rm n}$ ($n_{\rm p}$), its EOS $E(n,\delta)$ is usually expressed as the binding energy per nucleon as a function of the nucleon number density $n=n_{\rm n}+n_{\rm p}$ and the isospin asymmetry $\delta=(n_{\rm n}-n_{\rm p})/n$. The $E(n,\delta)$ can be expanded in terms of $\delta$ as

[TABLE]

where $E_{0}(n)=E(n,\delta=0)$ is the EOS of SNM, and the symmetry energy $E_{\rm sym}(n)$ is defined by

[TABLE]

It should be mentioned that the odd-order terms of $\delta$ vanish in Eq. (1) due to the exchange symmetry between protons and neutrons in nuclear matter. At the saturation density $n_{0}$, the $E_{0}(n)$ can be expanded in $\chi=(n-n_{0})/3n_{0}$ as

[TABLE]

where $E_{0}(n_{0})$ is the binding energy per nucleon of SNM at $n_{0}$, $K_{0}=\left.9n_{0}^{2}\frac{d^{2}E_{0}(n)}{dn^{2}}\right|_{n=n_{0}}$ is the incompressibility coefficient, and $J_{0}=\left.27n_{0}^{3}\frac{d^{3}E_{0}(n)}{dn^{3}}\right|_{n=n_{0}}$ is the skewness coefficient.

Around a reference density $n_{r}$, the $E_{\rm sym}(n)$ can be expanded in $\chi_{r}=(n-n_{r})/3n_{r}$ as

[TABLE]

where $L(n_{r})=\left.3n_{r}\frac{dE_{\rm sym}(n)}{dn}\right|_{n=n_{r}}$ is the density slope parameter and $K_{\rm sym}(n_{r})=\left.9n_{r}^{2}\frac{d^{2}E_{\rm sym}(n)}{dn^{2}}\right|_{n=n_{r}}$ is the density curvature parameter. At $n_{r}=n_{0}$, the $L(n_{r})$ and $K_{\rm sym}(n_{r})$ are reduced, respectively, to the well-known $L\equiv L(n_{0})$ and $K_{\rm sym}\equiv K_{\rm sym}(n_{0})$, which characterize the density dependence of the $E_{\rm sym}(n)$ around $n_{0}$.

2.2 The extended Skyrme-Hartree-Fock model

In this work, we use a single theoretical model, namely, the extended Skyrme-Hartree-Fock (eSHF) model (Chamel et al., 2009; Zhang & Chen, 2016) to simultaneously describe nuclear matter, finite nuclei and neutron stars. Compared to the standard SHF model (see, e.g., Ref. (Chabanat et al., 1997)), the eSHF model contains additional momentum- and density-dependent two-body forces to effectively simulate the momentum dependence of the three-body forces and can describe very well the properties of nuclear matter, finite nuclei and neutron stars (Zhang & Chen, 2016), which involve a wide density region from subsaturation to suprasaturation densities. We would like to emphasize that the density dependence of nuclear matter EOS and the $E_{\rm sym}(n)$ from the eSHF model is very flexible. In particular, within the eSHF model, the high-density $E_{\rm sym}(n)$ could be positive or negative while the $E_{\rm sym}(n)$ at saturation and subsaturation densities can be in nice agreement with the nuclear experimental constraints (Zhang & Chen, 2016). Accordingly, the eSHF model is especially suitable for our present motivation to explore the possibility for the existence of the supersoft high-density $E_{\rm sym}(n)$.

The extended Skyrme effective nucleon-nucleon interaction is taken to have a zero-range, density- and momentum-dependent form (Chamel et al., 2009; Zhang & Chen, 2016), i.e.,

[TABLE]

where we have $\bm{R}=(\bm{r}_{i}+\bm{r}_{j})/2$ and $\bm{r}=\bm{r}_{i}-\bm{r}_{j}$, $P_{\sigma}=(1+\bm{\sigma}_{i}\cdot\bm{\sigma}_{j})/2$ is the spin exchange operator, and $\bm{\sigma}_{i}$ ($\bm{\sigma}_{j}$) is the Pauli spin matrix. In addition, the relative momenta operators ${\bm{K}}=(\bm{\nabla}_{i}-\bm{\nabla}_{j})/2i$ and $\bm{K}^{\prime}=-(\bm{\nabla}_{i}-\bm{\nabla}_{j})/2i$ act on the right and left of the wave function, respectively. The interaction includes $14$ independent model parameters, i.e., the $13$ Skyrme force parameters $\alpha$, $t_{0}\sim t_{5}$, $x_{0}\sim x_{5}$, and the spin-orbit coupling constant $W_{0}$. The $13$ Skyrme force parameters can be expressed explicitly in terms of the following $13$ macroscopic quantities (pseudo-parameters) (Zhang & Chen, 2016): $n_{0}$, $E_{0}(n_{0})$, $K_{0}$, $J_{0}$, $E_{\rm sym}(n_{r})$, $L(n_{r})$, $K_{\rm sym}(n_{r})$, the isoscalar effective mass $m_{s,0}^{\ast}$, the isovector effective mass $m_{v,0}^{\ast}$, the gradient coefficient $G_{S}$, and the symmetry-gradient coefficient $G_{V}$, the cross gradient coefficient $G_{SV}$, and the Landau parameter $G_{0}^{\prime}$ of SNM in the spin-isospin channel. For the motivation of the present work, instead of directly using the $13$ Skyrme force parameters, it is very convenient to use the $13$ macroscopic quantities in the eSHF calculations for nuclear matter, finite nuclei and neutron stars, and the details can be found in Ref. (Zhang & Chen, 2016).

2.3 Tidal deformability of neutron stars

The tidal deformability (polarizability) $\lambda$ of NSs can be thought of as the NS fundamental f-modes with spherical harmonic index $l=2$ which can be treated as forced and damped harmonic oscillators driven by the external tidal field of the NS’s companion. The $\lambda$ is defined as the oscillation response coefficient (Flanagan & Hinderer, 2008), namely, the ratio of the neutron star’s quadrupole moment $Q_{ij}$ to the companion’s perturbing tidal field ${\cal E}_{ij}$ (in units with $c=G=1$ in this work) (Flanagan & Hinderer, 2008; Hinderer, 2008), i.e., $\lambda=-Q_{ij}/{\cal E}_{ij}$. The $\lambda$ is related to the dimensionless quadrupole tidal Love number $k_{2}$ and the NS radius $R$ by the relation $\lambda=\frac{2}{3}k_{2}R^{5}$. For a NS with mass $M$, the dimensionless tidal deformability $\Lambda_{M}$ is conventionally defined as

[TABLE]

The Love number $k_{2}$ depends on the details of the NS structure and it can be evaluated by (Hinderer, 2008)

[TABLE]

where $C=M/R$ is the NS compactness and $y=y(R)$ is determined by solving the following first-order differential equation:

[TABLE]

with

[TABLE]

In the above, $C_{s}^{2}\equiv dP(r)/d{\cal E}(r)$ is the squared sound speed of the NS matter. Eq. (8) for dimensionless $y(r)$ must be integrated with the general relativistic equations of hydrostatic equilibrium, namely, the famous Tolman-Oppenheimer-Volkoff (TOV) equations (Tolman, 1939; Oppenheimer & Volkoff, 1939):

[TABLE]

where $r$ is the radial coordinate, $M(r)$ is the enclosed mass inside the radius $r$, and ${\cal E}(r)$ ($P(r)$) is the energy density (pressure) at $r$. The boundary condition for Eq. (8) is $y(0)=2$ (Postnikov et al., 2010). For a given NS matter EOS $P({\cal E})$, one can calculate the NS mass $M$, radius $R$, Love number $k_{2}$, and $\Lambda_{M}$ with various central densities for the NS.

The NS contains core, inner crust, and outer crust. The density $n_{\rm{out}}$ separating the inner and outer crusts is taken to be $2.46\times 10^{-4}~{}\rm{fm}^{-3}$, and the core-crust transition density $n_{\rm t}$ is evaluated self-consistently by a dynamical approach (Xu et al., 2009). We assume here that the core is composed of $\beta$-stable and electrically neutral $npe\mu$ matter and its EOS can be calculated within the eSHF model. For the inner crust between densities $n_{\rm out}$ and $n_{\rm t}$, the EOS is constructed by interpolating with $P=a+b{\cal E}^{4/3}$ due to its complicated structure (Carrier et al., 2003). For the outer crust, we employ the well-known Baym-Pethick-Sutherland EOS in the density region of $6.93\times 10^{-13}~{}{\rm fm}^{-3}<n<n_{\rm out}$ and Feynman-Metropolis-Teller EOS for $n<6.93\times 10^{-13}~{}{\rm fm}^{-3}$ (Baym et al., 1971; Iida & Sato, 1997). The causality condition $dP/d{\cal E}\leq 1$ is guaranteed for all the NS calculations in the present work.

2.4 Fitting strategy for model parameters

In the eSHF model, there are totally $14$ model parameters, i.e., $n_{0}$, $E_{0}(n_{0})$, $K_{0}$, $J_{0}$, $E_{\rm sym}(n_{r})$, $L(n_{r})$, $K_{\rm sym}(n_{r})$, $m_{s,0}^{\ast}$, $m_{v,0}^{\ast}$, $G_{S}$, $G_{V}$, $G_{SV}$, $G_{0}^{\prime}$, and $W_{0}$. Following the same fitting strategy for model parameters as in Ref. (Zhou et al., 2019), we first fix $E_{\rm sym}(n_{c})=26.65$ MeV and $L(n_{c})=47.3$ MeV at the subsaturation density $n_{c}=0.11n_{0}/0.16$ according to the precise constraint $E_{\rm sym}(n_{c})=26.65\pm 0.2$ MeV (Zhang & Chen, 2013) by analyzing the binding energy difference of heavy isotope pairs and $L(n_{c})=47.3\pm 7.8$ MeV (Zhang & Chen, 2014) extracted from the electric dipole polarizability of 208Pb. In addition, the higher-order parameters $J_{0}$ and $K_{\rm sym}$ generally have small influence on the properties of finite nuclei but are crucial for the high-density nuclear matter EOS and the NS properties. To explore the effects of $J_{0}$ and $K_{\rm sym}$, we thus fix them at various values but with the other $10$ parameters being obtained by fitting the data of finite nuclei by minimizing the weighted sum of the squared deviations of the theoretical predictions from the experimental data, i.e.,

[TABLE]

where the $\bm{p}=(p_{1},...,p_{z})$ denote the $z$ dimensional model space, ${\cal O}_{i}^{({\rm th})}$ and ${\cal O}_{i}^{({\rm exp})}$ are the theoretical predictions and the corresponding experimental data, respectively, and $\Delta{\cal O}_{i}$ is the adopted error for balancing the relative weights of different types of observables (see, e.g., Ref. (Zhang & Chen, 2016)). The $30$ data of finite nuclei used in this work are listed in Table 1. As for $\Delta{\cal O}_{i}$, we use $1.0$ MeV and $0.01$ fm for the $E_{\rm{B}}$ and $r_{\rm{c}}$, respectively, and for the $E_{\rm{GMR}}$ we use the experimental error multiplied by $3.5$ to also consider the effect of the experimental error, while for the $\epsilon_{\rm{ls}}^{A}$ a $10\%$ relative error is employed. Considering the relatively larger uncertainty for $L(n_{c})=47.3\pm 7.8$ MeV (Zhang & Chen, 2014), we also investigate the cases with $L(n_{c})=39.5$ MeV and $55.1$ MeV.

3 Results and discussions

Using the fitting strategy described before, for $E_{\rm sym}(n_{c})=26.65$ MeV and $L(n_{c})=39.5$ MeV, $47.3$ MeV, and $55.1$ MeV, we construct a series of extended Skyrme parameter sets with fixed $J_{0}$ in the large range of ($-500,-300$) MeV and $K_{\rm sym}$ in ($-220,60$) MeV. As found in Ref. (Zhou et al., 2019), in order to be consistent with the constraint on the pressure of SNM in the density region of about $2n_{0}\sim 5n_{0}$ from the flow data in heavy-ion collisions (Danielewicz et al., 2002), the $J_{0}$ must be less than $-342$ MeV, i.e., the upper limit of $J_{0}$ is $J^{\rm up}_{0}=-342$ MeV, independent of the values of $L(n_{c})$ and $K_{\rm sym}$. Therefore, the flow data put strong constraint on the EOS of SNM at suprasaturation densities and can significantly limit the maximum mass of NSs (Zhou et al., 2019).

Shown in Fig. 1 is the NS maximum mass $M_{\rm max}$ vs $K_{\rm sym}$ using various extended Skyrme parameter sets. It is seen that for each $L(n_{c})$ with a fixed $J_{0}$, the $M_{\rm max}$ becomes insensitive to $K_{\rm sym}$ when the latter is larger than a critical value $K^{\rm crit}_{\rm sym}$. For $L(n_{c})=39.5$ MeV, $47.3$ MeV, and $55.1$ MeV, the value of $K^{\rm crit}_{\rm sym}$ is roughly $-130$ MeV, $-100$ MeV, and $-70$ MeV, respectively. These results imply that the $E_{\rm sym}(n)$ has little influence on the $M_{\rm max}$ when the $K_{\rm sym}$ is large enough. This can be understood from the fact that for the stiff high-density $E_{\rm sym}(n)$ with large $K_{\rm sym}$, the NS matter becomes almost isospin symmetric at high densities and the $M_{\rm max}$ hence essentially depends on the high-density EOS of SNM, which is mainly controlled by the $J_{0}$.

On the other hand, it is very interesting to see that for a fixed $J_{0}$, the $M_{\rm max}$ decreases drastically as the $K_{\rm sym}$ decreases when the $K_{\rm sym}$ is less than $K^{\rm crit}_{\rm sym}$. This means that the observed heaviest NS mass can rule out too soft high-density $E_{\rm sym}(n)$ with small $K_{\rm sym}$ values. From Fig. 1, one sees that for a fixed $K_{\rm sym}$, the $M_{\rm max}$ generally increases with $J_{0}$. Consequently, the extended Skyrme parameter sets with $J_{0}=J^{\rm up}_{0}=-342$ MeV generally predict the largest $M_{\rm max}$ in the eSHF model. For $L(n_{c})=(39.5,47.3,55.1)$ MeV, we obtain the largest $M_{\rm max}$ in the eSHF model as $(2.30,2.28,2.26)M_{\odot}$. Furthermore, we find for $L(n_{c})=(39.5,47.3,55.1)$ MeV, using the recently discovered heaviest NS with mass $2.14M_{\odot}$ sets a lower limit of $K_{\rm sym}$, namely, $K^{\rm low}_{\rm sym}=(-183,-157,-123)$ MeV, while using a NS maximum mass $2.01M_{\odot}$ gives $K^{\rm low}_{\rm sym}=(-198,-171,-134)$ MeV. Therefore, the existence of heavier NSs requires a stiffer lower bound of the high-density $E_{\rm sym}(n)$ with larger $K^{\rm low}_{\rm sym}$.

In addition, we note that for each $L(n_{c})$, the $\Lambda_{1.4}$ monotonically increases with $K_{\rm sym}$ ($J_{0}$) for a fixed $J_{0}$ ($K_{\rm sym}$) but the sensitivity on $K_{\rm sym}$ is much stronger than that on $J_{0}$ (Zhou et al., 2019). Therefore, the existence of the lower limit for $K_{\rm sym}$ (i.e., $K^{\rm low}_{\rm sym}$) leads to a lower bound of $\Lambda_{1.4}$, namely, $\Lambda^{\rm low}_{1.4}=(297,348,436)$ for $L(n_{c})=(39.5,47.3,55.1)$ MeV based on the so far measured heaviest NS mass $2.14M_{\odot}$. On the other hand, the lower bound of $\Lambda_{1.4}$ is found to be $\Lambda^{\rm low}_{1.4}=(245,286,391)$ for $L(n_{c})=(39.5,47.3,55.1)$ MeV by using a NS maximum mass $2.01M_{\odot}$. These results show that the lower bound of $\Lambda_{1.4}$ changes from $\Lambda^{\rm low}_{1.4}=286^{+105}_{-41}$ to $\Lambda^{\rm low}_{1.4}=348^{+88}_{-51}$ when the measured largest NS mass varies from $2.01M_{\odot}$ to $2.14M_{\odot}$. Therefore, the recently discovered heaviest NS, i.e., PSR J0740+6620 (Cromartie et al., 2019), puts a much stronger limit on $\Lambda^{\rm low}_{1.4}$, i.e., $\Lambda_{1.4}\geq 348^{+88}_{-51}$. The quite large lower bound of $\Lambda_{1.4}\geq 348^{+88}_{-51}$ combined with the upper limit $\Lambda_{1.4}\leq 580$ (Abbott et al., 2018) from the GW signal GW170817 leads to a stringent constraint on the $\Lambda_{1.4}$, i.e., $348^{+88}_{-51}\leq\Lambda_{1.4}\leq 580$. This will have important implications on the structure properties of NSs and the NS-involved GW detection in future.

Since the $\Lambda_{1.4}$ rapidly increases with $K_{\rm sym}$, the upper limit $\Lambda_{1.4}\leq 580$ from the GW signal GW170817 (Abbott et al., 2018) can set upper limits on $K_{\rm sym}$ for various values of $J_{0}$ as shown in Fig. 1. According to the allowed parameter space shown in Fig. 1, the recently discovered heaviest NS with mass $2.14M_{\odot}$ sets a upper limit of $K_{\rm sym}$, namely, $K^{\rm up}_{\rm sym}=(-46,-48,-53)$ MeV for $L(n_{c})=(39.5,47.3,55.1)$ MeV. We note that using a NS maximum mass $2.01M_{\odot}$ gives $K^{\rm up}_{\rm sym}=(-37,-39,-42)$ MeV for $L(n_{c})=(39.5,47.3,55.1)$ MeV. The existence of the upper and lower limits of $K_{\rm sym}$ can rule out too stiff and too soft high-density $E_{\rm sym}(n)$ and thus put strong constraints on the high-density behaviors of $E_{\rm sym}(n)$.

Figure 2 shows the density dependence of the symmetry energy according to the allowed parameter space for $J_{0}$ and $K_{\rm sym}$ with $L(n_{c})=(39.5,47.3,55.1)$ MeV as shown in Fig. 1. Fig. 2 (a) is obtained by using $2.01M_{\odot}$ as the NS maximum mass while Fig. 2 (b) is by using $2.14M_{\odot}$. Also included in Fig. 2 are the constraints on the $E_{\rm sym}(n)$ at subsaturation densities from midperipheral heavy-ion collisions of Sn isotopes (Tsang et al., 2009), the isobaric analog states (IAS) and combining the neutron skin data (IAS + NSkin) (Danielewicz & Lee, 2014), and the electric dipole polarizability ($\alpha_{D}$) in 208Pb (Zhang & Chen, 2015). For comparison, we further include in Fig. 2 (b) the results from some microscopic many-body approaches, namely, the non-relativistic Brueckner-Hartree-Fock (BHF) approach (Vidaña et al., 2009; Li et al., 2008), the relativistic Dirac-Brueckner-Hartree-Fock (DBHF) approach (Klähn et al., 2006; Sammarruca, 2010), and the variational many-body (VMB) approach (Akmal et al., 1998; Friedman & Pandaharipande, 1981; Wiringa et al., 1988). It is seen from Fig. 2 that the $E_{\rm sym}(n)$ with various values of $L(n_{c})$, $J_{0}$ and $K_{\rm sym}$ in the allowed parameter space are all in good agreement with the experimental constraints at subsaturation densities but exhibit very different high-density behaviors.

From Fig. 2 (a), one can see that in the case with a NS maximum mass $2.01M_{\odot}$, the lower bound of the $E_{\rm sym}(n)$ becomes negative when the density is larger than $n/n_{\rm nuc}\approx(5.6,6.3)$ for $L(n_{c})=(47.3,55.1)$ MeV (Here $n_{\rm nuc}=0.16$ fm*-3* represents nuclear normal density). We note the corresponding central density $n_{\rm cen}$ of the NS with mass $2.01M_{\odot}$ is $n_{\rm cen}/n_{\rm nuc}\approx(6.4,7.4)$ for $L(n_{c})=(47.3,55.1)$ MeV, indicating that the lower bound of the $E_{\rm sym}(n)$ already becomes negative at suprasaturation densities inside the NS, and therefore the corresponding $E_{\rm sym}(n)$ is supersoft, which can cause the appearance of a PNM core in the NS (Note: the higher-order symmetry energies, e.g., the fourth-order symmetry energy (Cai & Chen, 2012), may affect the proton fraction in NS matter, and especially in the case of the supersoft symmetry energy, they may obviously change the disappearance density of the proton fraction in NSs (Zhang & Chen, 2001)). Our results thus demonstrate that the supersoft high-density $E_{\rm sym}(n)$ can support a NS with mass $2.01M_{\odot}$, and at the same time it can describe very successfully the data of finite nuclei and agree well with the constraint from the flow data in heavy-ion collisions.

In the case with a NS maximum mass $2.14M_{\odot}$, on the other hand, it is very interesting to see from Fig. 2 (b) that the $E_{\rm sym}(n)$ is always positive and the supersoft $E_{\rm sym}(n)$ is clearly ruled out. This means that the eSHF model with a supersoft $E_{\rm sym}(n)$ cannot simultaneously describe the data of finite nuclei, the constraint on SNM EOS from flow data in heavy-ion collisions, and the NSs with mass $2.14M_{\odot}$. Our results therefore exclude the possibility for the appearance of a PNM core in NSs. Furthermore, while our results are consistent with most of the microscopic many-body calculations shown in Fig. 2 (b), they indeed rule out the VMB calculations with interactions WFF1 (i.e., AV14 plus UVII), WFF3 (i.e., UV14 plus TNI) (Wiringa et al., 1988) and FP (i.e., $v_{14}$ + TNI) (Friedman & Pandaharipande, 1981). Our present results also rule out many non-relativistic Skyrme and Gogny effective interactions that predict negative symmetry energy at suprasaturation densities (See, e.g., Refs. (Stone et al., 2003; Chen, 2017)). It is interesting to note that our results seem to support the relativistic mean-field description of nuclear matter, which generally cannot predict negative $E_{\rm sym}(n)$ at high densities due to the specific construction of meson exchanges (Chen et al., 2007; Dutra et al., 2014; Chen, 2017). Our present results on the constraints of high-density $E_{\rm sym}(n)$ may also have important implications on the poorly known effective three-body forces, short-range tensor forces and short-range nucleon-nucleon correlations (Xu & Li, 2010; Cai et al., 2018).

Finally, we would like to point out that including new degrees of freedom such as hyperons (Vidaña et al., 2011; Lonardoni et al., 2015), antikaon condensates (Gupta & Arumugam, 2013; Char & Banik, 2014), and quark matter (Bombaci et al., 2016; Alford & Sedrakian, 2017; Dexheimer et al., 2018) that could be present in the interior of NSs but neglected in the present work, usually softens the NS matter EOS, and in this case a stiffer high-density $E_{\rm sym}(n)$ would be necessary to obtain a NS with mass $2.14M_{\odot}$. Therefore, including the new degrees of freedom in NSs is also expected to rule out the supersoft high-density $E_{\rm sym}(n)$.

4 Conclusion

Within the theoretical framework of the eSHF model, we have demonstrated that a supersoft high-density symmetry energy cannot simultaneously describe the data of finite nuclei, the equation of state of symmetric nuclear matter at suprasaturation densities constrained from the flow data in heavy-ion collisions, and the maximum neutron star mass of $2.14M_{\odot}$, although it is still allowed if the maximum neutron star mass is $2.01M_{\odot}$. Therefore, the very recent discovery of PSR J0740+6620 rules out the supersoft high-density symmetry energy, which means it is unlikely to have a pure neutron matter core in neutron stars. Furthermore, we have found that the stiffer lower limit of the high-density symmetry energy based on the existence of $2.14M_{\odot}$ neutron stars leads to a quite large lower limit for $\Lambda_{1.4}$, i.e., $\Lambda_{1.4}\geq 348^{+88}_{-51}$, which is expected to have important implications on the future multimessenger observations of neutron-star-involved GW events.

The authors thank Tanja Hinderer, Bao-An Li and Zhen Zhang for useful discussions. This work was supported in part by the National Natural Science Foundation of China under Grant No. 11625521, the Major State Basic Research Development Program (973 Program) in China under Contract No. 2015CB856904, the Program for Professor of Special Appointment (Eastern Scholar) at Shanghai Institutions of Higher Learning, Key Laboratory for Particle Physics, Astrophysics and Cosmology, Ministry of Education, China, and the Science and Technology Commission of Shanghai Municipality (11DZ2260700).