# Impact of NICER’s Radius Measurement of PSR J0740+6620 on Nuclear Symmetry
Energy at Suprasaturation Densities

Nai-Bo Zhang11affiliation: School of Physics, Southeast University, Nanjing
211189, China and Bao-An Li22affiliation: Department of Physics and Astronomy,
Texas A$\&$M University-Commerce, Commerce, TX 75429, USA
∗Corresponding author<EMAIL_ADDRESS>∗

###### Abstract

By directly inverting several neutron star observables in the three-
dimensional parameter space for the Equation of State of super-dense neutron-
rich nuclear matter, we show that the lower radius limit for PSR J0740+6620 of
mass $2.08\pm 0.07~{}M_{\odot}$ from Neutron Star Interior Composition
Explorer (NICER)’s very recent observation sets a much tighter lower boundary
than previously known for nuclear symmetry energy in the density range of
$(1.0\sim 3.0)$ times the saturation density $\rho_{0}$ of nuclear matter. The
super-soft symmetry energy leading to the formation of proton polarons in this
density region of neutron stars is clearly disfavoured by the first radius
measurement for the most massive neutron star observed reliably so far.

Dense matter, equation of state, stars: neutron

##  1 Introduction

The PSR J0740+6620 having an updated mass of $2.08\pm 0.07~{}M_{\odot}$
(Fonseca et al., 2021) remains the most massive neutron star (NS) discovered
so far (Cromartie et al., 2019). Its radius of $13.7^{+2.6}_{-1.5}$ km (68%)
(Miller et al., 2021) or $12.39_{-0.98}^{+1.30}$ km (Riley et al., 2021) was
very recently inferred in two independent analyses of the X-ray data taken by
the Neutron Star Interior Composition Explorer (NICER) and the X-ray Multi-
Mirror (XMM-Newton) observatory. Combined with NICER’s earlier simultaneous
mass and radius measurement of PSR J0030+0451 (Miller et al., 2019; Riley et
al., 2019), the first radius measurement of the most massive NS has the strong
potential to reveal interesting new physics about the Equation of State (EOS)
of super-dense neutron-rich nuclear matter. Besides earlier predictions about
what uniquely new physics can be learned from the radii of massive NSs
compared to canonical ones, see, e.g., Xie & Li (2020); Han & Prakash (2020);
Drischler et al. (2021); Somasundaram & Margueron (2021), several new analyses
aiming at extracting new information about the EOS of super-dense matter from
NS observations including the latest NICER observations have already been
carried out (Biswas, 2021; Li et al., 2021; Raaijmakers et al., 2021; Pang et
al., 2021; Essick et al., 2021; Huth et al., 2021). Most of these studies are
based on the Bayesian statistical inference. Since the approximately 10-15%
statistical uncertainty on top of the roughly 10% systematic uncertainty (as
indicated by the difference of the results from the two analyses) of NICER’s
radius measurements are still quite large, improvements to the existing
constraints of nuclear EOS brought by the latest NICER measurement were found
to be generally small within the Bayesian statistical analyses.

After verifying the small effects of NICER’s new radius measurements on the
EOS in a Bayesian analysis (Xie & Li, 2021), here we perform a direct
inversion of the lower radius boundaries obtained by both Miller et al. (2021)
and Riley et al. (2021) in a three-dimensional (3D) high-density EOS parameter
space. Without the statistical averaging over the whole uncertainty range of
observables normally done in the Bayesian analyses, direct inversion of
specific values of NS observables can reveal interesting features that are
generally smeared out in the Bayesian statistical analyses. We skip the
inversion of the upper radius boundaries from the two NICER analyses as they
do not provide any tighter constraint on the EOS compared to what is available
from analyzing the upper limit of tidal deformation of GW170817. We show that
the 68% lower mass-radius boundaries from the two independent analyses by
Miller et al. (2021) and Riley et al. (2021) provide a much tighter lower
boundary than previously known for nuclear symmetry energy in the density
range of $(1.0\sim 3.0)\rho_{0}$, disfavoring the super-soft symmetry energy
necessary for the formation of proton polarons in this density region in NSs
and the associated phenomena predicted in the literature.

The rest of the paper is organized as follows. In the next section, we recall
the main challenges and ramifications of determining nuclear symmetry energy
at suprasaturation densities. We then summarize very briefly the direct
inversion approach of analyzing NS observables in the 3D high-density EOS
parameter space. We will demonstrate in Section 4 the scientific power of
NICER’s radius measurement of PSR J0740+6620 as the most massive NS with a
reliable mass in tightening the lower limit of high-density nuclear symmetry
energy. Effects of the remaining uncertainty of the slope parameter $L$ of
nuclear symmetry energy will be discussed in Section 5. Finally, we summarize
our main findings.

##  2 The challenges and ramifications of determining nuclear symmetry energy
at suprasaturation densities

The average energy per nucleon $E(\rho,\delta)$ in neutron-rich matter of
nucleon density $\rho=\rho_{n}+\rho_{p}$ and isospin asymmetry
$\delta\equiv(\rho_{n}-\rho_{p})/\rho$ can be written as (Bombaci & Lombardo,
1991)

$E(\rho,\delta)=E_{0}(\rho)+E_{\rm{sym}}(\rho)\cdot\delta^{2}+\mathcal{O}(\delta^{4})$
(1)

where $E_{0}(\rho)$ is the nucleon energy in symmetric nuclear matter (SNM)
while $E_{\rm{sym}}(\rho)$ is the nuclear symmetry energy. The latter measures
the energy cost to make nuclear matter more neutron-rich. It is essentially
the difference in energy per nucleon in pure neutron matter (PNM) and SNM. The
$E_{\rm{sym}}(\rho)$ is very poorly known especially at suprasaturation
densities mainly because of our poor knowledge about the spin-isospin
dependence of three-body nuclear forces, the isospin dependence of the tensor
force, and the related short-range nucleon-nucleon correlations in dense
neutron-rich matter (Li et al., 2018, 2019). It has been well known that
theoretical predictions of high-density $E_{\rm{sym}}(\rho)$ diverge broadly
at suprasaturation densities, ranging from large negative to positive values
at densities above about $(2\sim 3)\rho_{0}$ (Li et al., 1998; Steiner et al.,
2005; Baran et al., 2005; Li et al., 2008, 2014; Baldo & Burgio, 2016). While
there have been continuously a lot of new efforts and indeed some impressive
progresses in recent years in predicting nuclear symmetry energy at densities
below $2\rho_{0}$ using microscopic nuclear many-body theories, the situation
at higher densities remain essentially the same, see, e.g., the summary in
Figure 2 in the latest review (Burgio et al., 2021) or Figure 8 in (Li et al.,
2021). To our best knowledge, there is no fundamental physics principle
forbidding the $E_{\rm{sym}}(\rho)$ to become zero or even negative at
suprasaturation densities. In fact, based on variational nuclear many-body
theory calculations (Pandharipande & Garde, 1972; Wiringa et al., 1988), it
was pointed out already that when the repulsive short-range tensor force due
to the $\rho$ meson exchange in the isosinglet nucleon-nucleon interaction
channel in SNM becomes dominating at high densities, the energy in SNM
increases faster than that in PNM, leading to a decreasing
$E_{\rm{sym}}(\rho)$ with increasing density. The $E_{\rm{sym}}(\rho)$ may
then become zero or even negative above certain critical densities as the
$E_{\rm{sym}}(\rho)$ is approximately the energy difference between PNM (where
the tensor force in the isotriplet nucleon-nucleon interaction channel is
known to be negligible) and SNM. Indeed, it was demonstrated quantitatively by
Xu & Li (2010); Xu et al. (2013) that whether the $E_{\rm{sym}}(\rho)$
increases or decreases and when it may become negative depend strongly on
properties of the three-body force, the strength and short-range cut-off of
the tensor force as well as the related isospin dependence of short-range
nucleon-nucleon correlations.

The high-density behavior of $E_{\rm{sym}}(\rho)$ has significant
ramifications in both nuclear physics and astrophysics (Kutschera & Wójcik,
1993; Kutschera, 1994; Szmagliński et al., 2006; Baran et al., 2005; Li et
al., 2008; Ditoro et al., 2010; Baldo & Burgio, 2016; Zhou & Chen, 2019). For
example, variations of nuclear symmetry energy can lead to large changes in
the binding energy and spacetime curvature near the surface of NSs (Newton &
Li, 2009; He et al., 2015). To further test Einstein’s General Relativity in
the strong-field gravity regime against modified gravity theories or the
possible existence of a weakly interacting light boson mediating a new force
using massive NSs thus requires a reliable knowledge of the high-density
$E_{\rm{sym}}(\rho)$ (Krivoruchenko et al., 2009; Wen et al., 2009; Lin et
al., 2014; Jiang et al., 2015). Moreover, the density profile of proton
fraction in NSs at $\beta$-equilibrium is uniquely determined by the density
dependence of nuclear symmetry energy and has many consequences. For instance,
it determines the critical density above which the direct URCA process
responsible for the fast cooling of protoneutron stars can happen. In
particular, a super-soft $E_{\rm{sym}}(\rho)$, i.e., a decreasing/negative
symmetry energy at high densities may lead to some very interesting new
phenomena in the core of NSs. As the symmetry energy decreases with increasing
density, the proton fraction decreases correspondingly, leading to the
formation of proton polarons (isolated and localized single or small clusters
of protons surrounded by high-density neutrons) (Kutschera & Wójcik, 1993;
Kutschera, 1994) in the core of NSs. It can then strongly affect transport
coefficients of NS matter and can produce spontaneous magnetization in NSs
(Szmagliński et al., 2006). When the $E_{\rm{sym}}(\rho)$ becomes zero at very
high densities in the core of NSs, matter there are completely PNM. Since we
do not have a physical mechanism to handle negative symmetry energy in NSs in
the present work, we use a cut-off at zero symmetry energy at high densities
in cases where the $E_{\rm{sym}}(\rho)$ continuously decreases towards
negative values. Nevertheless, we notice that a negative symmetry energy and
its effects on NSs have been studies by others, see, e.g., Kutschera & Wójcik
(1993); Kutschera (1994); Szmagliński et al. (2006). In particular, it has
been shown that a negative $E_{\rm{sym}}(\rho)$ has dramatic effects on the
possible kaon condensation in the cores of NSs (Kubis & Kutschera, 2003).
Theoretically, in nuclear matter (not necessarily in NSs), if the
$E_{\rm{sym}}(\rho)$ becomes negative at super-high densities, the so-called
isospin separation instability happens (Kutschera, 1994; Li, 2002). Namely,
because of the $E_{\rm{sym}}(\rho)\cdot\delta^{2}$ term in the average nucleon
energy of Equation (1) in isospin-asymmetric nuclear matter, a uniform
symmetric nuclear matter ($\delta=0$) of energy $E_{0}(\rho)$ is unstable
energetically against being separated into regions of bulk pure neutron matter
($\delta=+1$) and pure proton matter ($\delta=-1$) when $E_{\rm{sym}}(\rho)$
is negative because then the
$E_{\rm{sym}}(\rho)(\delta=-1)^{2}+E_{\rm{sym}}(\rho)(\delta=+1)^{2}$ will be
negative instead of zero, making the total energy of the system less than the
$E_{0}(\rho)$ that a uniform SNM normally has.

To pin down the high-density behavior of $E_{\rm{sym}}(\rho)$ has been a major
scientific thrust of high-energy rare isotope beam facilities (U.S. LRP, 2015;
NuPECC LRP, 2017), see, e.g., Tsang et al. (2021); Russotto et al. (2021) for
examples of current and planned nuclear reaction experiments probing the high-
density $E_{\rm{sym}}(\rho)$. In fact, there were some circumstantial
evidences for a super-soft $E_{\rm{sym}}(\rho)$ at $\rho\geq 1.5\rho_{0}$ from
analyzing data of an earlier heavy-ion reaction experiment (Xiao et al.,
2009). However, the conclusion remains controversial (Feng, 2012a, b; Xie &
Zhang, 2014; Jhang et al., 2021) largely due to the strong model dependence
and the relatively small isospin asymmetry reached in heavy-ion reactions (Li,
2017). Thus, a tight lower boundary for the high-density $E_{\rm{sym}}(\rho)$
from measuring the radii of massive NSs are scientifically invaluable for both
nuclear physics and astrophysics.

While the radii and tidal deformabilities of canonical NSs provide some useful
constraints on the $E_{\rm{sym}}(\rho)$ around $(1\sim 2)\rho_{0}$, they do
not constrain significantly the $E_{\rm{sym}}(\rho)$ at densities above about
$2\rho_{0}$ as demonstrated clearly already in Zhang & Li (2019a); Xie & Li
(2019). This is mainly because the radii and tidal deformabilities of
canonical NSs are determined by the nuclear pressure around $(1\sim
2)\rho_{0}$ in these relatively light NSs (Lattimer & Prakash, 2001).
Moreover, the pressure in this density region is strongly controlled by the
density dependence of nuclear symmetry energy, see, e.g., Figure 5 in Li et
al. (2006) for a numerical demonstration. The radius measurement of canonical
NSs thus constrains mainly the symmetry energy in this density region. At
higher densities, the pressure is dominated by the SNM pressure. As pointed
out already by Lattimer & Prakash (2001), “there is a quantitative relation
between the radius and the pressure that does not depend upon the EOS at the
highest densities, which determines the overall softness or stiffness (and
hence, the maximum mass)”. Nowadays, when analyzing NS observables to extract
information about nuclear EOS, one always first applies the prerequisite that
all EOSs have to support the currently observed maximum mass of NSs. Under
this condition, the radius data of canonical NSs thus has little effect on the
SNM EOS in the whole density region but mainly constrain the
$E_{\rm{sym}}(\rho)$ below about $2\rho_{0}$. At densities above $2\rho_{0}$,
the remaining gap between the extracted upper and lower limits of
$E_{\rm{sym}}(\rho)$ remains very large, see, e.g., Figure 8 in the recent
review in Li et al. (2021), because the constrains on the $E_{\rm{sym}}(\rho)$
obtained from the radii of canonical NSs can only rule out some extreme cases
predicted by some theories. In fact, it was shown clearly before the GW170817
was discovered that to probe the high-density behavior of $E_{\rm{sym}}(\rho)$
one has to use the radii and/or tidal deformabilities of massive NSs while
those of the canonical ones only probe the $E_{\rm{sym}}(\rho)$ around
$\rho_{0}$, see, e.g., Fattoyev et al. (2013, 2014). NICER’s very recent
radius measurement of PSR J0740+6620 made this hope a reality for the first
time.

Figure 1: The symmetry energy (upper panels) and proton fraction (lower
panels) as functions of density from varying $L$ (left panel), $K_{\rm sym}$
(middle panel), and $J_{\rm sym}$ (right panel), respectively.

##  3 Direct inversion of neutron star observables in high-density EOS
parameter space

For completeness and ease of our discussions, here we summarize briefly the
main features of the NS EOS-metamodel we use in solving the NS inverse-
structure problem. More details can be found in our previous publications
(Zhang et al., 2018; Zhang & Li, 2019a, b, c, 2020; Xie & Li, 2019, 2020). We
assume the cores of NSs are made of totally charge neutral neutrons, protons,
electrons, and muons (the $npe\mu$ model) at $\beta$-equilibrium. Unlike the
widely used composition-degenerate spectral functions and/or piecewise
polytropes that directly parameterize the pressure as a function of energy or
baryon density, to probe the high-density symmetry energy we have to keep the
isospin dependence of the EOS and retain explicitly the composition
information at all densities. For this reason, we metamodel NS EOS by starting
at the single nucleon energy level with explicitly isospin dependence.
Specifically, we parameterize the $E_{0}(\rho)$ and $E_{\rm{sym}}(\rho)$
according to

$\displaystyle E_{0}(\rho)$ $\displaystyle=$ $\displaystyle
E_{0}(\rho_{0})+\frac{K_{0}}{2}(\frac{\rho-\rho_{0}}{3\rho_{0}})^{2}+\frac{J_{0}}{6}(\frac{\rho-\rho_{0}}{3\rho_{0}})^{3},$
(2) $\displaystyle E_{\rm{sym}}(\rho)$ $\displaystyle=$ $\displaystyle
E_{\rm{sym}}(\rho_{0})+L(\frac{\rho-\rho_{0}}{3\rho_{0}})$ (3)
$\displaystyle+$
$\displaystyle\frac{K_{\rm{sym}}}{2}(\frac{\rho-\rho_{0}}{3\rho_{0}})^{2}+\frac{J_{\rm{sym}}}{6}(\frac{\rho-\rho_{0}}{3\rho_{0}})^{3}$

where $E_{0}(\rho_{0})=-15.9\pm 0.4$ MeV is the binding energy and
$K_{0}\approx 230\pm 20$ MeV (Shlomo et al., 2006; Piekarewicz, 2010; Garg &
Colò, 2018) is the incompressibility at the saturation density $\rho_{0}$ of
SNM, while $E_{\rm sym}(\rho_{0})=31.7\pm 3.2$ MeV is the magnitude and
$L\approx 58.7\pm 28.1$ MeV is the slope of symmetry energy at $\rho_{0}$ (Li
& Han, 2013; Oertel et al., 2017) based on earlier surveys of over 50 analyses
of both terrestrial experiments and astrophysical observations, respectively.
A very recent survey of 24 new analyses of NS observations since GW179817
indicates that $L\approx 57.7\pm 19$ MeV and $K_{\rm{sym}}\approx-107\pm 88$
MeV at 68% confidence level (Li et al., 2021). In this study, we use $L=58.7$
as its most probable value and vary it within $\pm 20$ MeV. We keep the
$E_{0}(\rho_{0}),K_{0}$, and $E_{\rm sym}(\rho_{0})$ at their most probable
values given above as they have been relatively well determined and their
variations within their remaining uncertain ranges have been shown to have
little effects on the masses, radii, and tidal deformabilities of NSs Fattoyev
et al. (2013, 2014); Xie & Li (2019, 2020). As a reference, we notice that
very recent calculations based on a set of commonly used Hamiltonians
including two- and three-nucleon forces derived from chiral effective field
theory predicted that $L=47\pm 3$ MeV, $K_{\rm{sym}}=-146\pm 43$ MeV,
$J_{\rm{sym}}=90\pm 334$ MeV (Somasundaram et al., 2021). While in the same
calculations, with a prior $J_{0}=-300\pm 400$ MeV the posterior ranges
between $J_{0}=-573\pm 133$ MeV and $J_{0}=-172\pm 243$ MeV depending on the
scaling schemes used. The EOS parameter ranges we used are consistent with
these calculations.

The $K_{\rm{sym}}$, $J_{\rm{sym}}$, and $J_{0}$ are parameters characterizing
the EOS of super-dense neutron-rich nuclear matter. We directly invert NS
observables in the high-density EOS space allowed by the current uncertainties
of $K_{\rm{sym}}$, $J_{\rm{sym}}$, and $J_{0}$ as summarized in Zhang et al.
(2017). We note qualitatively that while the parameter $J_{0}$ (skewness of
SNM) controls the stiffness of SNM EOS at high densities, the $K_{\rm{sym}}$
(curvature of symmmetry energy) dominates the behavior of $E_{\rm sym}(\rho)$
around $2\rho_{0}$ and the parameter $J_{\rm{sym}}$ (skewness of symmetry
energy) controls the $E_{\rm sym}(\rho)$ at densities above $(2\sim
3)\rho_{0}$. The slope parameter $L$ dominates the behavior of $E_{\rm
sym}(\rho)$ around $\rho_{0}$ but does not have much effects on its high-
density behavior. It is also known from previous Bayesian analyses and direct
inversions that the $K_{\rm{sym}}$ and $J_{\rm{sym}}$ affect most strongly on
the radii of massive NSs while the $L$ parameter influences most strongly the
radii of canonical NSs (Zhang & Li, 2019b; Xie & Li, 2020).

The roles of the parameters $L$, $K_{\rm{sym}}$, and $J_{\rm{sym}}$ mentioned
above are illustrated in Figure 1 where the $E_{\rm sym}(\rho)$ (upper panels)
and the corresponding proton fraction (lower panels) in NSs at $\beta-$
equilibrium are shown as functions of density by varying the slope $L$ (left
panel), curvature $K_{\rm sym}$ (middle panel), and skewness $J_{\rm sym}$
(right panel), respectively. They can sample essentially all currently known
variations of density dependence of nuclear symmetry energy. We loosely refer
the $E_{\rm sym}(\rho)$ that quickly decreases to zero above $(2-3)\rho_{0}$
as “super-soft”. This terminology has been used by the nuclear physics
community, see, e.g., Xiao et al. (2009). The $E_{\rm sym}(\rho)$ that
increases much faster than linear is loosely referred as “super-stiff”, based
on the fact that the original Relativistic Mean Field (RMF) model predicts a
linear increase with density due to the $\rho-$meson exchange and the RMF
prediction is often referred as being “stiff”. As shown in Figure 1,
combinations of large negative (positive) $K_{\rm{sym}}$ and/or $J_{\rm{sym}}$
values lead to super-soft (super-stiff) symmetry energies at high densities.
For a super-soft $E_{\rm{sym}}(\rho)$, the proton fraction $X_{p}$ quickly
decreases to zero, indicating the formation of PNM. Conversely, the $X_{p}$
quickly approaches 1/2 indicative of a SNM at super-high densities for the
super-stiff $E_{\rm{sym}}(\rho)$.

It is worth emphasizing that the Equations (2) and (3) for $E_{0}(\rho)$ and
$E_{\rm{sym}}(\rho)$ have dual meanings that can be confusing. If the
functions $E_{0}(\rho)$ and $E_{\rm{sym}}(\rho)$ are known apriori as in all
forward-modelings and in some inversions using energy density density
functionals, the Equations (2) and (3) are Taylor expansions up to the third
order in $\chi\equiv(\rho-\rho_{0})/3\rho_{0}$. It is well known that the
Taylor expansion may not converge at high densities. In fact, it has been
shown very recently in Cai & Li (2021) that even if all the coefficients in
the Equations (2) and (3) are well determined by whatever experiments and/or
theories, the Taylor expansions even up to $\chi^{5}$ can not reproduce
accurately the known $E_{0}(\rho)$ and $E_{\rm{sym}}(\rho)$ at densities
around $(3-4)\rho_{0}$. To remedy this problem, an auxiliary function approach
instead of the original Taylor expansion was proposed. However, it has not
been implemented in any analysis of NS properties yet.

On the other hand, in our metamodeling of NS EOS in this work and our previous
works using both Bayesian and direct inversion approaches, the $E_{0}(\rho)$
and $E_{\rm{sym}}(\rho)$ are NOT known apriori. The Equations (2) and (3) for
$E_{0}(\rho)$ and $E_{\rm{sym}}(\rho)$ are simply parameterizations with their
coefficients to be inferred directly from data. They are thus not Taylor
expansions of any known functions. But they are purposely being parameterized
in forms like Taylor expansions so that we can use existing knowledge about
the coefficients of Taylor expansions of nuclear EOS as we mentioned above as
guidances in setting the range and selecting the prior probability
distribution functions (PDFs) of the parameters involved. Mathematically, the
parameterizations approach asymptotically to Taylor expansions at the limit of
$\rho\rightarrow\rho_{0}$. Indeed, simply referring the parameters $L$,
$K_{\rm sym}$, and $J_{\rm sym}$ as the slope, curvature, and skewness of
symmetry energy is misleading. But neglecting this subtle difference has been
a convenient and common practice in the literature.

Similar to the question of how many segments of polytropes one should use in
parameterizing the EOS at high densities to achieve model independent results,
to what order of $\chi$ one can practically use and how the results may be
different are interesting questions. Some answers to these questions were
given within the Bayesian approach in Xie & Li (2019) and Xie & Li (2021a).
For example, the kurtosis parameter $Z_{0}$ is the coefficient of the
$\chi^{4}$ term in Taylor expanding the SNM EOS $E_{0}(\rho)$. Its theoretical
value is known roughly to be between -5330 MeV and +5038 MeV (Margueron et
al., 2018; Antić et al., 2019). Using the empirical pressure in SNM from
nuclear collective flow in relativistic heavy-ion collisions as data
(Danielewicz et al., 2002), the most probable values and 68% confidence
boundaries of both the inferred $K_{0}$ and $J_{0}$ vary appreciable
especially for the latter depending on whether or not the $Z_{0}$ term is
included in parameterizing the $E_{0}(\rho)$ (Xie & Li, 2021a). While the
Bayesian inferred $Z_{0}$ value of $600\pm 200$ MeV is much more accurate than
its prior value from the literature, its uncertainty remains very large and it
affects mostly the SNM EOS above about $(3-4)\rho_{0}$. As another example, as
demonstrated in Section 4 of Xie & Li (2019), the most probable values and
confidence boundaries of $L$ and $K_{\rm sym}$ extracted from NS observables
depend significantly on whether the $J_{\rm{sym}}$ and $J_{0}$ terms are
considered in the EOS parameterizations. Indeed, in the analyses by different
groups in the literature, often the $J_{\rm{sym}}$ is neglected as its range
is very poorly known as we mentioned earlier. In fact, whether the
$J_{\rm{sym}}$ is considered or not is partially responsible for the
variations of $L$ and $K_{\rm sym}$ reported in the literature as pointed out
already in Xie & Li (2019). As to the kurtosis parameter $Z_{\rm{sym}}$ of
symmetry energy, it is currently known to be around $-1800\pm 800$ MeV
(Somasundaram et al., 2021) and is not considered in our work.

For visualization purposes, our direct inversion as we shall illustrate next
is limited to 3D. Given the above information and for the purpose of
extracting the $E_{\rm{sym}}(\rho)$ mostly in the region of $(2-4)\rho_{0}$,
we perform the direct inversion in the $J_{0}-K_{\rm sym}-J_{\rm sym}$ 3D
high-density EOS parameter space. For this purpose, the parameterizations of
$E_{0}(\rho)$ and $E_{\rm{sym}}(\rho)$ are sufficient. It is not meaningful
and impractical to include even higher order terms before the $J_{\rm sym}$
parameter is better determined. Using less terms is also not meaningful as the
$J_{\rm sym}$ is most important for determining the $E_{\rm{sym}}(\rho)$
around $(2-4)\rho_{0}$. Thus, within the direct inversion approach we
currently can not answer the question about how our results may change if more
or less terms are used in parameterizing the EOS, and as we discussed earlier
there is currently no theoretical consensus on the $E_{\rm{sym}}(\rho)$ above
$2\rho_{0}$ for us to compare with. Therefore, all of our results should be
understood within the context and cautions discussed above.

The pressure in the $npe\mu$ core of NSs is calculated from

$P(\rho,\delta)=\rho^{2}\frac{d\epsilon(\rho,\delta)/\rho}{d\rho},$ (4)

where
$\epsilon(\rho,\delta)=\epsilon_{n}(\rho,\delta)+\epsilon_{l}(\rho,\delta)$
denotes the total energy density with contributions from both nucleons
($\epsilon_{n}(\rho,\delta)$) and leptons ($\epsilon_{l}(\rho,\delta)$). To
see clearly where the symmetry energy comes in, we note in particular that the
nucleon energy density $\epsilon(\rho,\delta)=\rho E(\rho,\delta)+\rho M$ with
$M$ being the average mass of nucleons. Given the density dependence of
nuclear symmetry energy $E_{\rm sym}(\rho)$, the density profile of isospin
asymmetry $\delta(\rho)$ (or the corresponding proton fraction $x_{p}(\rho)$)
is determined by the $\beta$-equilibrium and charge neutrality conditions
$\mu_{n}-\mu_{p}=\mu_{e}=\mu_{\mu},~{}\rho_{p}=\rho_{e}+\rho_{\mu}$ in terms
of the chemical potential
$\mu_{i}=\frac{\partial\epsilon(\rho,\delta)}{\partial\rho_{i}}$ and particle
density $\rho_{i}$ of particle $i$.

We connect the core EOS described above with the NV EOS (Negele & Vautherin,
1973) for the inner crust and the BPS EoS (Baym et al., 1971) for the outer
crust. The crust-core transition density and pressure are determined
consistently for each EOS generated for the core as discussed in detail in
Zhang et al. (2018).

Figure 2: Constant surfaces of NS observables and causality condition in the
3-dimensional $K_{\rm sym}-J_{\rm sym}-J_{0}$ EOS parameter space when $L$ is
fixed at its currently known most probable value of $L=58.7$ MeV: the NS
minimum maximum mass of $M=2.01$ M⊙ (green surface), the dimensionless tidal
deformability of canonical NS: $\Lambda_{1.4}=580$ at 90% confidence level
(yellow surface) and $\Lambda_{1.4}=427$ at 68.3% confidence level (pink
surface), and the causality surface (blue) on which the sound speed equals the
speed of light in centers of the most massive NSs. The red arrows point to the
corresponding surfaces while the black arrows show directions that satisfy the
corresponding constraints.

Figure 3: The central density (left panel) and radius of a NS with the maximum
mass of 2.01 M⊙ as functions of $K_{\rm sym}$ and $J_{\rm sym}$ parameters of
the symmetry energy.

The $npe\mu$ model outlined above can be considered as the minimum model for
NSs. Other particles and new phases of matter may appear especially in massive
NSs. However, there are many uncertainties and unknowns associated with these
new particles and phases. In particular, to our best knowledge, there is
currently no community consensus regarding the existence and type of a hadron-
quark phase transition but a hadron-quark duality in understanding properties
of NSs. For example, strong evidences for a hadron-quark phase transition were
found for massive but not canonical NSs in Annala et al. (2020). However, a
very recent analysis including NICER’s latest observation for PSR J0740+6620
disfavours the presence of a strong first-order phase transition from nuclear
matter to exotic forms of matter, such as quark matter, inside NSs (Pang et
al., 2021).

The Bayesian analysis of canonical NS data in Xie & Li (2021b) using the above
metamodel EOS for hadronic matter and the constant speed of sound (CSS) model
for quark matter (Alford et al., 2013) strongly favors a large quark core even
in canonical NSs with the most probable hadron-quark transition density at
$\rho_{t}/\rho_{0}=1.6^{+1.2}_{-0.4}$ at 68% confidence level. Similarly, low
hadron-quark transition densities were also found very recently in two
independent Bayesian analyses using similar NS data (Tang et al., 2021; Li et
al., 2021). While the PDFs of some of the six hadronic EOS parameters in
Equations (2) and (3) we extracted with and without considering the first-
order hadron-quark phase transition are mildly different, they are all
compatible with current constraints within the uncertain ranges of these
parameters discussed earlier. In particular, as shown in Figure 12 in Li et
al. (2021), while the 68% confidence boundaries of $E_{\rm{sym}}(\rho)$ at
suprasaturation densities become wider when the phase transition is considered
as three more parameters in the CSS model are included in the Bayesian
analyses of the same set of NS data, existing constrains from nuclear
theories, terrestrial experiments, and astrophysical observations can not
distinguish the results obtained with and without considering the hadron-quark
phase transition. Thus, there is clearly a hadron-quark duality in describing
NS observables similar to the situation in relativistic heavy-ion collisions.
Namely, purely hadronic models and the ones considering a hadron-quark phase
transition can describe equally well all observations within currently
acceptable/existing uncertainties of models/observations. Since we are not
trying to find ways to break this duality, the $npe\mu$ model is sufficient
for the purposes of this work. In our opinion, the choice of using the
$npe\mu$ model is reasonable given the current situation of the field although
it explores a limited EOS parameter space for NSs.

Unlike in Bayesian statistical analyses normally using the Markov-Chain Monte
Carlo sampling, our direct inversions of NS observables in the 3D high-density
EOS parameter space are done by brute force. Namely, using the metamodel EOS
outlined above, for a given observable of NSs we loop through the
3-dimensional high-density EOS parameter space spanned by the $K_{\rm{sym}}$,
$J_{\rm{sym}}$, and $J_{0}$ parameters within their currently known
uncertainties given earlier. At each step of the three loops, we call the
Tolman-Oppenheimer-Volkov (TOV) solver. In looping through the entire high-
density EOS parameter space, we also require all accepted EOSs to satisfy the
causality condition, retain the dynamical stability through out the entire NS
as well as a positive pressure at the crust-core transition density. All
accepted EOSs giving the same NS observable will be presented by a constant
observable surface in the $K_{\rm{sym}}$-$J_{\rm{sym}}$-$J_{0}$ EOS parameter
space. We emphasize that every point on the constant surface represents a
unique EOS and all EOSs on the surface have the same probability to reproduce
the specified NS observable. Crosslines of two observable surfaces and/or
those between an observable surface and the causality surface will set
boundaries for the acceptable high-density EOS parameter space.

As an example, shown in Figure 2 are the direct inversions of (1) the upper
limit of the tidal deformability of NSs in GW170817 (Abbott et al., 2018)
$\Lambda_{1.4}=580$ at 90% confidence level (yellow surface) and
$\Lambda_{1.4}=427$ at 68% confidence level (pink surface), respectively, (2)
the causality surface (blue surface) on which the speed of sound equals the
speed of light ($v^{2}_{s}=dP/d\epsilon=c^{2}$) at the central density of the
most massive NS supported by the nuclear pressure at each point with the
specific EOS there (Zhang & Li, 2019a), and (3) the constant maximum mass of
$M$=2.01 M⊙ which is the lower limit of the newly updated mass of PSR
J0740+6620 (Fonseca et al., 2021; Miller et al., 2021; Riley et al., 2021). In
this work, we use this mass as the minimum maximum mass of NSs that all EOSs
have to be able to support. Namely, all EOSs have to generate stable mass-
radius curves having maximum masses higher than or equal to 2.01 M⊙. These
constant surfaces were found by calling the TOV solver at each step in looping
through the $K_{\rm{sym}}$-$J_{\rm{sym}}$-$J_{0}$ EOS parameter space in the
ranges specified by the three axises of the plot.

The $M=2.01$ M⊙ (green) surface together with the causality surface (blue)
limit the allowed range of the $J_{0}$, thus the high-density SNM EOS. They
are essentially vertical to the $J_{0}$ axis when the $E_{\rm{sym}}(\rho)$ is
stiff with large positive $K_{\rm{sym}}$ and $J_{\rm{sym}}$ values. This is
because when the symmetry energy is stiff, as shown in Figure 1, the isospin
asymmetry $\delta=1-2X_{p}$ at $\beta$-equilibrium is small (Quantitatively
from solving the equation $\mu_{n}-\mu_{p}=\mu_{e}=\mu_{\mu}\approx 4\delta
E_{\rm{sym}}(\rho)$. Qualitatively, it can be understood simply from
minimizing the $E_{\rm{sym}}(\rho)\cdot\delta^{2}$ term in Equation (1)). The
nuclear pressure is then dominated by the contribution from SNM EOS. However,
when the $E_{\rm{sym}}(\rho)$ becomes super-soft with large negative
$K_{\rm{sym}}$ and/or $J_{\rm{sym}}$ values, nuclear matter becomes very
neutron-rich with large $\delta$ values approaching 1. Then, the symmetry
energy contribution to nuclear pressure becomes large and can be even
negative. Thus to support the same NS mass $M=2.01$ M⊙, the SNM EOS has to
become stiffer with higher $J_{0}$ values, leading to the bending up of the
$M=2.01$ M⊙ surface at the right-back corner where the $K_{\rm{sym}}$ and/or
$J_{\rm{sym}}$ are largely negative. To see the maximum density reached in
this NS of mass $M$=2.01 M⊙ and its radius, shown in Figure 3 are its central
density (left panel) and radius (right panel) as functions of $K_{\rm sym}$
and $J_{\rm sym}$ parameters of the symmetry energy. The central density
ranges between about $(5.3\sim 9.5)\rho_{0}$ while the corresponding radius is
between about $(8.5\sim 12.3)$ km. As we shall discuss next, the absolutely
maximum mass of non-rotating NSs is determined by the causality surface. If a
NS with mass $M=2.01$ M⊙ is below the maximum mass supported by a given EOS,
its central density will be lower and its radius larger than the constant
surfaces shown in Figure 3 as indicated by the two black arrows.

For most of the $J_{\rm{sym}}$ values the causality surface bends downward as
the $K_{\rm{sym}}$ parameter goes from positive to large negative values. As
discussed in detail in Zhang & Li (2019a), this is because both the NS maximum
mass $M_{\rm{cau}}$ on the causality surface that the EOS can support and the
corresponding radius $R_{\rm{cau}}$ decrease when the $E_{\rm{sym}}(\rho)$
becomes more soft as the $K_{\rm{sym}}$ decreases. Since the maximum density
$\rho_{\rm{max}}$ reached in this maximum mass NS scales with
$M_{\rm{cau}}/R^{3}_{\rm{cau}}$, it increases very quickly with the decreasing
curvature $K_{\rm{sym}}$ towards super-soft $E_{\rm{sym}}(\rho)$. Therefore,
the causality condition $v_{s}=c$ can be reached at much smaller $J_{0}$
values as the $K_{\rm{sym}}$ decreases. We notice that the variation of the
causality surface with the skewness parameter $J_{\rm{sym}}$ is more
complicated as its effects on both the maximum mass and the corresponding
radius are non-monotonic depending on the SNM EOS parameter $J_{0}$.

Continuing the discussions of the results shown in Figure 2, the crossline
between the causality surface and the $\Lambda_{1.4}=580$ (427) surface
determines the upper limit of $E_{\rm{sym}}(\rho)$ at 90% (68%) confidence
level. Since the crosslines between the causality surface and the upper radius
limits from NICER’s observation of PSR J0740+6620 are on the left of the
$\Lambda_{1.4}=580$ (90% confidence) surface, they are not shown here. We
emphasize that before considering the lower radius limits from NICER’s
observation of PSR J0740+6620, the crossline between the causality surface and
the NS minimum maximum mass $M=2.01$ M⊙ sets the lower boundary for the high-
density symmetry energy. Both the 68% upper and lower boundaries of
$E_{\rm{sym}}(\rho)$ from these crosslines will be compared to the new 68%
lower boundaries set by the lower radius boundaries from the two independent
analyses of NICER’s observation of PSR J0740+6620 in the next section.

##  4 The power of radius measurement of massive neutron stars on constraining
high-density nuclear EOS

Shown in Figure 4 are the inversions of the NICER’s radius lower limit
$R_{1.28}=11.52$ km (pink surface) for PSR J0030+0451 (Riley et al., 2019) as
well as $R_{2.01}=11.41$ (orange surface) km (Riley et al., 2021) and
$R_{2.01}=12.2$ km (yellow surface) (Miller et al., 2021) for PSR J0740+6620
together with the causality and the minimum maximum mass $M=2.01$ M⊙ surfaces
in the 3D high-density EOS parameter space. While the $R_{1.28}=11.52$ km
surface is on the right of the $R_{2.01}=11.41$ km surface, it is shown here
to make a comparison. Similarly, while the $M$=2.01 M⊙ surface is below the
$R_{2.01}=11.41$ km surface and thus does not provide the lower boundary of
$E_{\rm sym}(\rho)$ any more, it is shown here to compare the constraint on
$E_{\rm sym}(\rho)$ obtained when only the mass of PSR J0740+6620 is measured
with that when both its mass and radius are measured. The black arrows
indicate the direction to the allowed EOS space.

Figure 4:  Similar to Figure 2 but replacing the tidal deformations with
NICER’s radius lower limit $R_{1.28}=11.52$ km (pink surface) for PSR
J0030+0451 (Riley et al., 2019) as well as $R_{2.01}=11.41$ (orange surface)
km (Riley et al., 2021) and $R_{2.01}=12.2$ km (yellow surface) (Miller et
al., 2021) for PSR J0740+6620.

Several interesting and important observations can be made from comparing the
constant surfaces of NS observables and the causality condition as well as
their crosslines shown in Figure 4. In particular, we emphasize the following
points:

  * •

The power in limiting the EOS parameter space using simultaneously the
knowledge about both the mass and radius of a massive NS can be seen clearly
from comparing the $M=2.01$ M⊙ surface and the $R_{2.01}=12.2$ or 11.41 km
surface. The separation between the $R_{2.01}=12.2$ km and $R_{2.01}=11.41$ km
surfaces indicate the importance of further reducing the uncertainty of the
radius measurement. For the discussions here we use the $R_{2.01}=12.2$ km
surface, while the lower boundaries of $E_{\rm{sym}}(\rho)$ using both radii
will be obtained and compared later. There is a huge gap in the direction of
$J_{0}$ between the two surfaces in the area where the $E_{\rm{sym}}(\rho)$ is
super-soft. The gap closes gradually as the $E_{\rm{sym}}(\rho)$ becomes
super-stiff towards the front-left corner. This is understood again because
the pressure is dominated by the SNM EOS when the $E_{\rm{sym}}(\rho)$ is
super-stiff and the corresponding isospin asymmetry $\delta$ vanishes. In this
region of the high-density EOS parameter space, the only mechanism to change
both the NS mass and its radius simultaneously is through the variation of
$J_{0}$.

On the other hand, it is seen that in the space where the $E_{\rm{sym}}(\rho)$
is super-soft, the $R_{2.01}=12.2$ km surface is almost vertical while the
$M=2.01$ M⊙ surface is still rather flat (which can be better seen in Figure
2), meaning that the variation of $J_{0}$ has essentially no effect on the
radius of PSR J0740+6620 in this EOS parameter region. Most importantly, the
crossline between the $R_{2.01}=12.2$ km surface and the causality surface
sets the new lower boundary of $E_{\rm{sym}}(\rho)$. Compared to the lower
boundary set by the crossline between the $M=2.01$ M⊙ and causality surfaces
shown in Figure 2, it moved upward (become stiffer) significantly as we shall
discuss more quantitatively below. Therefore, the most important effect of
knowing simultaneously both the radius and mass of the most massive NS
observed is in pinning down the high-density symmetry energy especially its
lower boundary.

  * •

The power in limiting the EOS parameter space of knowing the radius of a
massive NS compared to that of a light one can be clearly seen by comparing
the $R_{2.01}=12.2$ km surface and the $R_{1.28}=11.52$ km surface. The latter
is rather vertical because in the whole EOS space considered, all values of
$J_{0}$ can support a NS of mass $M=1.28$ M⊙ that has a relatively low central
density and the $J_{0}$ has almost no effect on the radii of light NSs, in
contrast to the situation for the most massive NS observed so far as we
discussed above. The separation between these two surfaces clearly indicates
the strong sensitivity of NS radius to the variation of $E_{\rm{sym}}(\rho)$.
Since both surfaces are the lower radius limits of the two NSs considered, it
is obvious that the $R_{2.01}=12.2$ km surface provides a more tight lower
limit for the $E_{\rm{sym}}(\rho)$, while the crossline of the
$R_{1.28}=11.52$ km surface with either the causality or the
$M=2.01~{}M_{\odot}$ surface gives a lower boundary very close to that from
the crossline between the causality and $M=2.01$ $M_{\odot}$ surfaces.

Figure 5: The projection of the crossline between surfaces of causality
condition and $M_{\rm max}=2.01$ M⊙ (green line), $R_{2.01}=11.41$ km (red
line) and $R_{2.01}=12.2$ km (blue line), respectively, on the $K_{\rm
sym}-J_{\rm sym}$ plane for $L=58.7$ MeV. The arrows point to the directions
satisfying the corresponding constraints.

Figure 6: Upper panel: The constrained lower limit of symmetry energy as a
function of density for $L=58.7$ MeV when only mass is observed (green line)
or both mass and radius are observed (red and blue lines). The constrained
upper limit from $\Lambda_{1.4}=427$ is shown as orange line. Lower panel: The
corresponding lower limit of proton fraction as a function of density. The
horizontal line shows the critical proton fraction 11.1% for the direct URCA
process to occur. The arrows point to the directions satisfying the
corresponding constraints.

We now turn to the 68% confidence boundaries of the high-density symmetry
energy from the crosslines discussed above. The function $E_{\rm{sym}}(\rho)$
is determined by the $K_{\rm sym}$ and $J_{\rm sym}$ parameters when the $L$
and $E_{\rm sym}(\rho_{0})$ are fixed at their most probable values given
earlier. Shown in Figure 5 are projections of the crossline between surfaces
of causality condition and $M_{\rm max}=2.01$ M⊙ (green line),
$R_{2.01}=11.41$ km (red line), and $R_{2.01}=12.2$ km (blue line),
respectively, on the $K_{\rm sym}-J_{\rm sym}$ plane for $L=58.7$ MeV. As
discussed above, these are the lower boundaries. The arrows point to the
directions satisfying the corresponding constraints. We note that along these
boundaries the values of $J_{0}$ are different as shown in Figures 2 and 4. To
see again effects of simultaneously measuring the mass and radius of massive
NSs, it is interesting to note that the boundary from the crossline between
the causality and $M=2.01~{}M_{\odot}$ surfaces sets the lower limit of
$K_{\rm sym}$ at -349 MeV. The information about the radius lower limit moves
the latter significantly. For example, the lower radius limit $R_{2.01}=12.2$
km moves it to about -215 MeV, thus reducing the uncertainty range of $K_{\rm
sym}$ by a factor of about 1.6. However, it is seen that the difference due to
the systematic error in analyzing NICER’s radius measurements is quite large.
Moreover, the $J_{\rm sym}$ parameter is still not constrained by the NS
observations available. Hopefully it can be constrained by future radius
measurements of even more massive NSs, messengers directly from the core of
NSs, e.g., neutrinos, postmerger high-frequency gravitational waves and/or
properties of hyper/super-massive remnants following mergers of two NSs.

Using the Equation (3), the 68% confidence boundaries in the $K_{\rm
sym}-J_{\rm sym}$ plane can be translated into constraints for the symmetry
energy $E_{\rm sym}(\rho)$ as shown in the upper panel of Figure 6. For a
comparison, the upper 68% confidence boundary from the crossline between the
$\Lambda_{1.4}=427$ and causality surfaces is shown as the orange line. In
translating the boundaries in the $K_{\rm sym}-J_{\rm sym}$ parameter plane to
the lower/upper boundaries of $E_{\rm{sym}}(\rho)$, since each point on the
boundary generates a new function $E_{\rm{sym}}(\rho)$ we have to compare all
functions generated along a given boundary to find the lowest/highest
$E_{\rm{sym}}(\rho)$ value at each density $\rho$. Then, linking all
lowest/highest values at all densities gives us the lower/upper boundary of
symmetry energy as a function of density. Because of the finite bin size in
density used in the translation, the resulting $E_{\rm{sym}}(\rho)$ functions
are not very smooth at all densities. Clearly, the measured radius of the most
massive NS observed so far improves significantly our knowledge about nuclear
symmetry energy especially at densities higher than about $2.5\rho_{0}$. More
quantitatively, the lowest symmetry energies at $2\rho_{0}$ and $3\rho_{0}$
and the corresponding proton fractions from the three crosslines are listed in
Table 1. Again, the lower radius limits from analyzing NICER’s radius
measurement of PSR J0740+6620 significantly increases the allowed lowest value
of symmetry energy at these two densities.

As mentioned earlier, the density profile of proton fraction $X_{P}(\rho)$
determined uniquely by the $E_{\rm sym}(\rho)$ is a critical quantity for
addressing several interesting NS physics issues. The constrained proton
fraction $X_{P}(\rho)$ is shown in the lower panel of Figure 6. The direct
URCA process is predicted to occur once the proton fraction exceeds the
critical fraction of about $11.1\%$ (shown as the horizontal dashed line)
(Klähn et al., 2006). The results of our analyses indicate that the direct
URCA process is disfavoured below about $3\rho_{0}$ but remains uncertain at
higher densities.

Table 1: The 68% confidence lower limits of symmetry energy and proton fraction at 2$\rho_{0}$ and 3$\rho_{0}$ from the crosslines of the causality condition and the three NS observables indicated in the first column. Causality & | $E_{\rm sym}(2\rho_{0})$ | $E_{\rm sym}(3\rho_{0})$ | $X_{p}(2\rho_{0})$ | $X_{p}(3\rho_{0})$
---|---|---|---|---
| (MeV) | (MeV) | (%) | (%)
$M_{\rm max}=2.01$ M⊙ | 37.32 | 20.96 | 4.1 | 0.5
$R_{2.01}=11.41$ km | 38.16 | 38.46 | 4.6 | 3.5
$R_{2.01}=12.2$ km | 44.26 | 50.96 | 6.7 | 7.2

Moreover, the possible formation of proton polarons and several related
phenomena as well as the trigger of isospin separation instability in the core
of NSs all depend critically on the proton fraction at high
densities(Kutschera & Wójcik, 1993; Kutschera, 1994; Kubis & Kutschera, 2003;
Szmagliński et al., 2006). One common condition is that the symmetry energy
has to start decreasing with increasing density above certain critical
density, leading to the gradual disappearance of protons in the core of NSs.
From Figure 6, it is seen that the lowest boundary (the green line determined
by the crossline between the causality and the $M=2.01~{}M_{\odot}$ surfaces)
of $X_{p}$ permits this at densities higher than about $2\rho_{0}$ when only
the mass $M=2.01~{}M_{\odot}$ is observed for PSR J0740+6620. In fact, many
phenomenological and microscopic models predicted high-density behaviors of
the $E_{\rm sym}$ and the corresponding $X_{p}$ similar to the dashed line,
see, e.g., the 10 realistic EOSs used in Szmagliński et al. (2006) in studying
the proton polaron formations in NSs. Interestingly, however, the lower
boundary of $X_{p}$ from using the lower radius limit $R_{2.01}=12.2$ km or
11.41 km disfavours the formation of proton polarons at least upto about
$3.0\rho_{0}$. While it is still more likely at higher densities, our present
model assuming NSs are made of $npe\mu$ matter may break down and thus can not
be used reliably to make predictions at higher densities.

Figure 7: Comparisons of the constant surfaces of $R_{1.28}=11.52$ km (upper
panel) and $R_{2.01}=12.2$ km (lower panel) for $L=40$, $58.7$, and $80$ MeV,
respectively.

##  5 Effects of the remaining uncertainty of the slope parameter $L$ of
nuclear symmetry energy

Figure 8: Comparisons of the constant surfaces of mass $M=2.01$ M⊙ (upper
panel) and causality condition (lower panel) for $L=40$, $58.7$, and $80$ MeV,
respectively.

Figure 9: The combined constraints on the lower limit of symmetry energy for
$L$ = 40 (dash lines), 58.7 (solid lines), and 80 (dash-dot-dot lines) MeV.

In the above discussions, the slope $L$ of symmetry energy is set at its
fiducial value of $L=58.7$ MeV, while it is known that it still has an
uncertainty of about 20 MeV based on the latest survey of available analyses
in the literature (Li et al., 2021). The $L$ parameter characterizes mainly
the symmetry energy around $\rho_{0}$. It is well known that $L$ has strong
imprints on the radii and tidal deformations of canonical NSs. However, it is
unclear how it may affect the radii and masses of massive NSs as well as the
causality surface. Most importantly, how does $L$ affect the lower boundaries
of the symmetry energy we extracted above? To address these questions we have
re-calculated everything with $L=40$ and 80 MeV and compared the results with
the ones we obtained above with $L=58.7$ MeV. As examples, shown in Figure 7
are comparisons of the constant surfaces of $R_{1.28}=11.52$ km (upper panel),
and $R_{2.01}=12.2$ km (lower panel) with $L=40$, $58.7$, and $80$ MeV,
respectively. As expected, the increase of $L$ shifts the surface
$R_{1.28}=11.52$ km to smaller $K_{\rm sym}$ values. For the light NSs, when
$L$ becomes larger, $K_{\rm sym}$ has to become smaller to reproduce the same
radii. This is consistent with our previous findings shown in Zhang & Li
(2019b). Interestingly, however, the above expectation is also true for the
radii of massive NSs only in the space where the symmetry energy is super-soft
as indicated by the surfaces of $R_{2.01}=12.2$ km for different $L$ in the
lower panel. In the region where the symmetry energy is super-stiff, the $L$
is seen to have little effect on the $R_{2.01}=12.2$ km surface. As we
discussed earlier, when the symmetry energy is stiff (larger $K_{\rm sym}$
and/or $J_{\rm sym}$), the pressure is dominated by the SNM EOS, namely the
$J_{0}$ parameter. Then both the mass and radius are mainly determined by the
$J_{0}$ with little influence from the symmetry energy. Moreover, since the
$K_{\rm sym}$ and $J_{\rm sym}$ are both so big in the super-stiff region, the
$L$ parameter thus has little effect on the $R_{2.01}=12.2$ km surface when
the symmetry energy is super-stiff. For similar reasons, as shown in Figure 8,
the variation of $L$ has little effect on the constant surface of $M=2.01$ M⊙
and the causality surface, except in the space where the symmetry energy is
super-soft.

Since the crosslines determining the boundaries of $E_{\rm sym}$ are not in
the super-soft region, effects of $L$ on these boundaries are expected to be
small or moderate. Summarized in Figure 9 are constraints on the 68%
confidence lower boundaries of symmetry energy from the two analyses of the
NICER data using $L=40$, $58.7$, and $80$ MeV, respectively. It is seen that
below about 3$\rho_{0}$ the variation of the lower boundary of $E_{\rm sym}$
due to variation of $L$ by $\pm 20$ MeV is compatible with the variation due
to the systematic error of NICER’s radius analysis.

##  6 Summary and conclusions

In summary, by directly inverting several NS observables in the high-density
EOS parameter space we have shown that the lower radius limits from NICER’s
very recent observation of PSR J0740+6620 set a much tighter lower boundary
than previously known for nuclear symmetry energy in the density range of
$\sim(1.0-3.0)\rho_{0}$. The super-soft symmetry energy leading to the
formation of proton polarons and the related phenomena predicted to occur in
this density region in the core of NSs is clearly disfavoured by NICER’s
radius measurement for the most massive NS observed so far.

Our work has two major caveats. Firstly, we did not consider any kind of phase
transition although the central density may reach far above the normally
expected hadron-quark transition density in PSR J0740+6620. Naturally, the
additional flexibility from adding a phase transition widens the range of
hadronic matter EOS parameters that agrees with the lower radius bound of PSR
J0740+6620. Moreover, the physical meaning of nuclear symmetry energy ends at
the hadron-quark phase transition. Secondly, our EOS parameterizations ends at
the $\chi^{3}$ terms. Including higher-order terms adds more flexibilities of
describing the high-density behavior of nuclear EOS. Since our direct
inversions of NS observables by brute force are limited to 3D, effects of
high-order terms are not investigated in this work. These may introduce some
systematic uncertainties in our results. Our findings should thus be
understood with these cautions in mind. Nevertheless, the relative positions
of the lower boundaries of high-density symmetry energy from using different
NS observables extracted from the same approach are physically sound and
qualitatively reliable.

## Acknowledgement

We would like to thank Wen-Jie Xie for helpful discussions. This work is
supported in part by the U.S. Department of Energy, Office of Science, under
Award Number DE-SC0013702, the CUSTIPEN (China-U.S. Theory Institute for
Physics with Exotic Nuclei) under the US Department of Energy Grant No. DE-
SC0009971, the National Natural Science Foundation of China under Grant No.
12005118, and the Shandong Provincial Natural Science Foundation under Grants
No. ZR2020QA085.

## References

  * Abbott et al. (2018) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2018, Phys. Rev. Lett., 121, 161101
  * Alford et al. (2013) Alford, M.G., Han, S., & Prakash, M. 2013, Phys. Rev. D, 88, 083013
  * Antić et al. (2019) Antić, S., Chatterjee, D., Carreau, T., et al. 2019, JPhG, 46
  * Annala et al. (2020) Annala, E., Gorda, T., Kurkela, A. et al. 2020, Nat. Phys., 16, 907
  * Baldo & Burgio (2016) Baldo, M., & Burgio, G. F. 2016, Prog. Part. Nucl. Phys., 91, 203
  * Baran et al. (2005) Baran, V., Colonna, M., Greco, V., & Di Toro, M. 2005, Phys. Rep., 410, 335
  * Baym et al. (1971) Baym, G., Pethick, C. J., & Sutherland, P. 1971, ApJ, 170, 299
  * Biswas (2021) Biswas, B. 2021, arXiv:2105.02886
  * Bombaci & Lombardo (1991) Bombaci, I., & Lombardo, U. 1991, Phys. Rev. C, 44, 1892
  * Burgio et al. (2021) Burgio, G. F., Vidana, I., Schulze, H.-J., et al. 2021, arXiv:2105.03747, Progress in Particle and Nuclear Physics in press.
  * Cai & Li (2021) Cai, B. J. & Li, B. A. 2021, Phys. Rev. C, 103, 054611.
  * Cromartie et al. (2019) Cromartie, H. T., Fonseca, E., Ransom, S. M., et al. 2019, NatAs, 4, 72
  * Danielewicz et al. (2002) Danielewicz, P., Lacey, R., & Lynch, W. G. 2002, Science, 298, 1592
  * Ditoro et al. (2010) Ditoro, M., Baran, V., Colonna, M., & Greco, V. 2010, J. Phys. G., 37, 083101
  * Drischler et al. (2021) Drischler, C., Han, S., Lattimer, J. M., et al. 2021, Phys. Rev. C, 103, 045808
  * Essick et al. (2021) Essick, R., Landry, P., Schwenk, A., et al. 2021, arXiv:2107.05528
  * Fattoyev et al. (2013) Fattoyev, F. J., Carvajal, J., Newton, W. G., & Li, B. A. 2013, Phys. Rev. C, 87, 015806
  * Fattoyev et al. (2014) Fattoyev, F. J., Newton, W. G., & Li, B. A. 2014, Phys. Rev. C, 90, 022801
  * Feng (2012a) Feng, Z. Q. 2012a, Nucl. Phys. A, 878, 3
  * Feng (2012b) Feng, Z. Q. 2012b, Physics Letters B 707, 83
  * Fonseca et al. (2021) Fonseca, E., Cromartie, H. T., Pennucci, T. T., et al. 2021, arXiv:2104.00880
  * Garg & Colò (2018) Garg, U., & Colò, G. 2018, Prog. Part. Nucl. Phys., 101, 55
  * Han & Prakash (2020) Han, S. & Prakash, M. 2020, ApJ, 899, 164
  * He et al. (2015) He, X. T., Fattoyev, F. J., Li, B. A., & Newton, W. G. 2015, Phys. Rev. C, 91, 015810
  * Huth et al. (2021) Huth, S., Pang, P. T. H., Tews, I., et al. 2021, arXiv:2107.06229
  * Jhang et al. (2021) Jhang, G., Estee, J., Barne, J., et al. 2021, Physics Letters B, 813, 136016
  * Jiang et al. (2015) Jiang, W. Z., Li, B. A., & Fattoyev, F. J. 2015, EPJA, 51, 119
  * Krivoruchenko et al. (2009) Krivoruchenko, M. I., Simkovic F., & Faessler, A. 2009, Phys. Rev. D, 79, 125023
  * Klähn et al. (2006) Klähn, T., Blaschke, D., Typel, S., et al. 2006, Phys. Rev. C, 74, 035802
  * Kubis & Kutschera (2003) Kubis, S., & Kutschera, M. 2003, Nucl. Phys. A, 720, 189
  * Kutschera & Wójcik (1993) Kutschera, M., & Wójcik, W. 1993 Phys. Rev. C, 47, 1077
  * Kutschera (1994) Kutschera, M. 1994, Physics Letters B, 340, 1
  * Lattimer & Prakash (2000) Lattimer, J. M., & Prakash, M. 2000, Phys. Rep., 333, 121
  * Lattimer & Prakash (2001) Lattimer, J. M., & Prakash, M. 2001, ApJ, 550, 426
  * Li et al. (2021) Li, A., Miao, Z., Han, S., & Zhang, B. 2021, ApJ, 913, 27
  * Li et al. (1998) Li, B. A., Ko, C. M., & Bauer W. 1998, Int. Jour. Mod. Phys. E, 7, 147
  * Li (2002) Li, B. A. 2002, Phys. Rev. Lett., 88, 192701
  * Li et al. (2006) Li, B.A., Chen, L.W., Ko, C.M., & Steiner, A.W, 2006, Revista mexicana de f sica, 52, 56 http://www.scielo.org.mx/pdf/rmf/v52s4/v52s4a11.pdf
  * Li et al. (2008) Li, B. A., Chen, L. W., & Ko, C. M. 2008, Phys. Rep., 464, 113
  * Li & Han (2013) Li, B. A., & Han, X. 2013, Phys. Lett. B, 727, 276
  * Li et al. (2014) Li, B. A., Ramos, À., Verde G., & Vidaña, I. (Eds.) 2014, EPJA, 50, 9
  * Li (2017) Li, B. A. 2017, Nuclear Physics News, 27, 7
  * Li et al. (2018) Li, B. A., Cai, B. J., Chen, L.W., & Xu J. 2018, Progress in Particle and Nuclear Physics, 99, 29
  * Li et al. (2019) Li, B. A., Krastev, P. G., Wen, D. H., & Zhang, N. B. 2019, EPJA, 55, 117
  * Li et al. (2021) Li, B. A., Cai, B. J., Xie, W. J., & Zhang, N. B. 2021, Universe 7, 182
  * Lin et al. (2014) Lin, W., Li, B. A., Chen, L. W., Wen, D. H., & Xu, J. 2014, JPhG, 41, 075203
  * Margueron et al. (2018) Margueron, J., Hoffmann Casali, R., & Gulminelli, F. 2018, Phys. Rev. C, 97, 025805
  * Miller et al. (2019) Miller, M. C., Lamb, F. K., Dittmann, A. J., et al. 2019, ApJ, 887, L24
  * Miller et al. (2021) Miller, M. C., Lamb, F. K., Dittmann, A. J., et al. 2021, arXiv:2105.06979
  * Negele & Vautherin (1973) Negele, J. W., & Vautherin, D. 1973, Nucl. Phys. A, 207, 298
  * Newton & Li (2009) Newton, W. G., & Li, B. A. 2009, Phys. Rev. C, 80, 065809
  * Oertel et al. (2017) Oertel, M., Hempel, M., Klähn, T., & Typel, S. 2017, Rev. Mod. Phys., 89, 015007
  * Pandharipande & Garde (1972) Pandharipande, V. R., & Garde, V. K. 1972, Physics Letters B, 39, 608
  * Pang et al. (2021) Pang, P. T. H., Tews, I., Coughlin, M. W., et al. 2021, arXiv:2105.08688
  * Piekarewicz (2010) Piekarewicz, J. 2010, JPhG, 37, 064038
  * Raaijmakers et al. (2021) Raaijmakers, G., Greif, S. K., Hebeler, K., et al. 2021, arXiv:2105.06981
  * U.S. LRP (2015) Reaching for the Horizon, the 2015 Long Range Plan for Nuclear Science, https://www.osti.gov/biblio/1296778-reaching-horizon-long-range-plan-nuclear-science
  * Riley et al. (2019) Riley, T. E., Watts, A. L., Bogdanov, S., et al. 2019, ApJ, 887, L21
  * Riley et al. (2021) Riley, T. E., Watts, A. L., Ray, P. S., et al. 2021, arXiv:2105.06980
  * Russotto et al. (2021) Russotto, P., Le Fèvre, A., Łukasik, J., et al. 2021, arXiv:2105.09233
  * Shlomo et al. (2006) Shlomo, S., Kolomietz, V. M., & Colò G. 2006, EPJA, 30, 23
  * Somasundaram et al. (2021) Somasundaram, R., Drischler, C., Tews, I., et al. 2021, Phys. Rev. C, 103, 045803.
  * Somasundaram & Margueron (2021) Somasundaram, R., & Margueron, J. 2021, arXiv:2104.13612
  * Steiner et al. (2005) Steiner, A. W., Prakash, M., Lattimer, J. M., & Ellis, P. J. 2005, Phys. Rep., 410, 325
  * Szmagliński et al. (2006) Szmagliński, A., Wójcik, W., & Kutschera, M. 2006, Acta Phys. Pol. B 37, 277
  * NuPECC LRP (2017) The Nuclear Physics European Collaboration Committee (NuPECC), Long Range Plan 2017 Perspectives in Nuclear Physics, https://www.esf.org/fileadmin/user_upload/esf/Nupecc-LRP2017.pdf
  * Tang et al. (2021) Tang, S.P., Jiang, J.L., Gao, W.H., Fan, Y.Z., & Wei, D.M. 2021, Phys. Rev. D, 103, 063026
  * Tsang et al. (2021) Tsang, M. B., et al. 2021, The Symmetry Energy Project, https://groups.nscl.msu.edu/hira/sepweb/pages/about.html
  * Wen et al. (2009) Wen, D. H., Li, B. A., & Chen, L. W. 2009, Phys. Rev. Lett., 103, 211102
  * Wiringa et al. (1988) Wiringa, R. B., Fiks, V., & Fabrocini, A. 1988, Phys. Rev. C, 38, 1010
  * Worley et al. (2008) Worley, A., Krastev, P. G., & Li, B. A. 2008, ApJ, 685, 390
  * Xiao et al. (2009) Xiao, Z. G., Li, B. A., Chen, L. W., Yong, G. C., & Zhang, M. 2009, Phys. Rev. Lett., 102, 062502
  * Xie & Zhang (2014) Xie, W. J., & Zhang, F. S. 2014, Physics Letters B, 735, 250
  * Xie & Li (2019) Xie, W. J., & Li, B. A. 2019, ApJ, 883, 174
  * Xie & Li (2020) Xie, W. J., & Li, B. A. 2020, ApJ, 899, 4
  * Xie & Li (2021a) Xie, W. J., & Li, B. A. 2021, J. of Phys. G 48, 025110
  * Xie & Li (2021b) Xie, W.J., & Li, B.A. 2021, Phys. Rev. C, 103, 035802
  * Xie & Li (2021) Xie, W. J., & Li, B. A. 2021 in preparation.
  * Xu & Li (2010) Xu, C. & Li, B. A. 2010, Phys. Rev. C, 81, 064612
  * Xu et al. (2013) Xu, C., Li, A., & Li, B. A. 2013, Journal of Physics: Conference Series, 420, 012090
  * Zhang et al. (2017) Zhang, N. B., Cai, B. J., Li, B. A., Newton, W. G., & Xu, J. 2017, Nucl. Sci. Tech., 28, 181
  * Zhang et al. (2018) Zhang, N. B., Li, B. A., & Xu, J. 2018, ApJ, 859, 90
  * Zhang & Li (2019a) Zhang, N. B., & Li, B. A. 2019, EPJA, 55, 39
  * Zhang & Li (2019b) Zhang, N. B., & Li, B. A. 2019, JPhG, 46, 014002
  * Zhang & Li (2019c) Zhang, N. B., & Li, B. A. 2019, ApJ, 879, 99
  * Zhang & Li (2020) Zhang, N. B., & Li, B. A. 2020, ApJ, 883, 61
  * Zhou & Chen (2019) Zhou, Y., & Chen, L.W., 2019, ApJ, 886, 52