# Central Speed of Sound, Trace Anomaly and Observables of Neutron Stars from
Perturbative Analyses of Scaled TOV Equations

Bao-Jun<EMAIL_ADDRESS>Quantum Machine Learning Laboratory, Shadow
Creator Inc., Shanghai 201208, China  Bao-An<EMAIL_ADDRESS>Department of Physics and Astronomy, Texas A$\&$M University-Commerce,
Commerce, TX 75429-3011, USA  Zhen<EMAIL_ADDRESS>Sino-
French Institute of Nuclear Engineering and Technology, Sun Yat-Sen
University, Zhuhai 519082, China

###### Abstract

The central speed of sound (SS) measures the stiffness of the Equation of
State (EOS) of superdense neutron star (NS) matter. Its variations with
density and radial coordinate in NSs in conventional analyses often suffer
from uncertainties of the specific nuclear EOSs used. Using the central SS and
NS mass/radius scaling obtained from solving perturbatively the scaled Tolman-
Oppenheimer-Volkoff (TOV) equations, we study the variations of SS, trace
anomaly and several closely related properties of NSs in an EOS-model
independent manner. We find that the SS increases with the reduced central
pressure $\widehat{P}_{\textrm{c}}\equiv
P_{\textrm{c}}/\varepsilon_{\textrm{c}}$ (scaled by the central energy density
$\varepsilon_{\textrm{c}}$), and the conformal bound for SS tends to break
down for NSs with masses higher than about 1.9$M_{\odot}$. The ratio
$P/\varepsilon$ is upper bounded as $P/\varepsilon\lesssim 0.374$ around the
centers of stable NSs. We demonstrate that it is an intrinsic property of
strong-field gravity and is more relevant than the perturbative QCD bound on
it. While a sharp phase transition at high densities characterized by a sudden
vanishing of SS in cores of massive NSs are basically excluded, the
probability for a continuous crossover signaled by a peaked radial profile of
SS is found to be enhanced as $\widehat{P}_{\textrm{c}}$ decreases, implying
it likely happens near the centers of massive NSs. Moreover, a new and more
stringent causality boundary as $R_{\max}/\textrm{km}\gtrsim
4.73M_{\textrm{NS}}^{\max}/M_{\odot}+1.14$ for NS M-R curve is found to be
excellently consistent with observational data on NS masses and radii.
Furthermore, new constraints on the ultimate energy density and pressure
allowed in NSs before collapsing into black holes are obtained and compared
with earlier predictions in the literature.

###### pacs:

21.65.-f, 21.30.Fe, 24.10.Jv

##  I Introduction and Conclusions

The speed of sound squared (SSS) is defined as
$s^{2}=\mathrm{d}P/\mathrm{d}\varepsilon$ [1] where $P$ and $\varepsilon$ are
respectively the pressure and energy density of the matter under
consideration. The conditions of stability and causality together require that
$0\leq s^{2}\leq 1$ (adopting $c=1$). Due to the superdense nature of neutron
star (NS) matter [2, 3, 4, 5, 6, 9, 10, 12, 11, 7, 8, 13, 14, 15, 16, 17, 18,
19, 20, 21, 22, 23, 24, 25] especially near its core [26, 27, 28, 29, 30, 31,
32, 33, 34, 35, 36], the speed of sound (SS) can become very close to 1 [37]
or even exceed it in certain models of dense matter. As a measure of the
stiffness of nuclear EOS, the behavior of SS is closely related to possible
phase transitions (PTs) and/or a continuous crossover in NSs. In particular, a
sharp vanishing of SSS (i.e., $s^{2}=0$ or a constant $P$ in a finite range of
$\varepsilon$) signals the occurrence of a first-order PT while a smooth
reduction of $s^{2}$ indicates a continuous crossover. Several possible
mechanisms for the PTs and/or crossover in NSs have been proposed in the
literature, see, e.g., Refs. [38, 39, 40, 41, 42, 43, 44, 45] for recent
reviews. While possible imprint of the SS on NS observables has been studied
extensively, the SS in NSs is not a quantity directly measurable and there are
still many interesting issues about it to be addressed.

Much progress has been made in understanding the internal structures and
observational properties of NSs during the last decade thanks especially to
the new opportunities provided by the era of multimessenger astronomy since
the discovery of GW170817 [46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58,
61, 64, 59, 60, 65, 66, 67, 62, 68, 63]. For instance, the NICER (Neutron Star
Interior Composition Explorer) and XMM-Newton simultaneously measured the
mass-radius for PSR J0740+6620 [69, 70, 71, 72] and that for PSR J0030+0451
[73, 74]. However, the relevant phases of matter in NS cores and therefore the
EOSs there still remain ambiguous [75, 76, 77, 78, 79, 80, 84, 81, 82, 83, 85,
86, 87, 95, 94, 88, 89, 90, 91, 92, 93, 96, 97, 98, 99, 100, 101, 102, 103,
104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114]. This is partially
because of the lack of direct accesses to the core EOS without relying on the
still uncertain structures and EOSs in other parts of NSs especially the NS
crust. Thus, many interesting and critical questions regarding properties of
superdense NS matter remain to be addressed. To outline our motivations and
set the context of this work, we list below a few examples:

  1. (a)

What is the EOS of densest visible matter existing in the Universe? Can it be
accessed/constrained directly using certain astrophysical data such as
observed NS radii and/or masses without using any EOS model?

  2. (b)

There exists an ultimate limit for the pressure/energy density in NSs above
which the system becomes unstable. What is this ultimate limit? Can the
observational data on NSs tell this quantity?

  3. (c)

Is the limit on the ratio $P/\varepsilon\leq 1$ (from the principle of
causality) relevant/sufficient in NS cores? If it is insufficient, how and to
which extend can this limit be improved? In fact, $P\leq\varepsilon$ is
equivalent to $s^{2}\leq 1$ only for a linear EOS $P\propto\varepsilon$, while
the EOS in NSs (especially in their cores) can be significantly nonlinear. As
a result, the upper bound for $P/\varepsilon$ in NSs is expected to be smaller
than 1. What are the implications of this refined upper bound for
$P/\varepsilon$ on certain NS properties which are determined by the balancing
acts of gravitational attraction and nuclear repulsion? Is this refined limit
a fundamental bound from General Relativity (GR) theory of gravity or
intrinsically related to microscopic theories of dense matter EOS such as
perturbative QCD (pQCD)?

  4. (d)

How can one understand the possible violation of the conformal bound (CB,
i.e., SSS is less than 1/3 based on the pQCD prediction) for SSS in massive
NSs? Can one calculate/estimate certain characteristic quantities related to
the conformality via the general-relativistic NS structural (TOV) equations
alone without using any EOS model for NS matter?

  5. (e)

Is the available causality boundary in the literature for NSs on the mass-
radius (M-R) curve a tight one? Can one refine or improve it? (Such boundary
should be consistent with NS observational data besides microscopic physics
constraints). Similarly, how large can the compactness mass/radius of a NS be?

  6. (f)

What is the radial profile of SS in NSs? Can one locate its peak if it exists?

Recently, we studied the question (a) above using a novel approach [115] which
facilitates extracting the central EOS of NS matter in the maximum-mass
configuration (where
$\mathrm{d}M_{\textrm{NS}}/\mathrm{d}\varepsilon_{\textrm{c}}=0$ with
$\varepsilon_{\textrm{c}}$ being the central energy density on the mass-radius
curve) directly from the observational data without relying on any specific
EOS model. The method was developed from analyzing perturbatively the
dimensionless Tolman-Oppenheimer-Volkoff (TOV) equations for NS internal
variables scaled by $\varepsilon_{\textrm{c}}$. In particular, a formula for
SSS in the center of these NSs was obtained as,

$s_{\textrm{c}}^{2}=\widehat{P}_{\textrm{c}}\left(1+\frac{1}{3}\frac{1+3\widehat{P}_{\textrm{c}}^{2}+4\widehat{P}_{\textrm{c}}}{1-3\widehat{P}_{\textrm{c}}^{2}}\right),$
(1)

here $\widehat{P}_{\textrm{c}}\equiv P_{\textrm{c}}/\varepsilon_{\textrm{c}}$
is the ratio of the central pressure $P_{\textrm{c}}$ over the central energy
density $\varepsilon_{\textrm{c}}$. From Eq. (1), one immediately finds that
$s_{\textrm{c}}^{2}>0$, therefore excluding a sharp PT (at high densities)
occurring in cores of maximum-mass configuration NSs. Whether it allows a
continuous crossover in the cores depends on the behavior of the SSS away from
the center as sketched in panel (a) of FIG. 1. A continuous crossover
indicates $s^{2}(\widehat{\varepsilon})>s_{\textrm{c}}^{2}$ (panel (b) where
the SSS at the center is smaller than that at a finite distance away) and a
possible peak in $s^{2}$ somewhere away from the center as the SSS should
become zero at the core-crust transition density, here $\widehat{r}$ is the
dimensionless distance from the center (Section II),
$\widehat{\varepsilon}=\varepsilon/\varepsilon_{\textrm{c}}$ and $D$ is a
coefficient characterizing the variation of $s^{2}$ with $\widehat{r}$
(Section VII). It is seen from Eq. (1) that the chance for
$s^{2}(\widehat{\varepsilon})>s_{\textrm{c}}^{2}$ is reduced as
$\widehat{P}_{\textrm{c}}$ increases since $s_{\textrm{c}}^{2}$ is upper
bounded. Therefore, for NSs with their $\widehat{P}_{\textrm{c}}$’s close to
its upper bound of about 0.374 from causality [115], the crossover near the
center is unfavored with high probabilities. In the limiting case of
$s_{\textrm{c}}^{2}\to 1$, the probability naturally approaches zero and it is
probably in this situation $s_{\textrm{c}}^{2}>s^{2}(\widehat{\varepsilon})$
(panel (c) of FIG. 1). Moreover, we can deduce from Eq. (1) that
$s_{\textrm{c}}^{2}\to 1/3$ occurs earlier than $\widehat{P}_{\textrm{c}}\to
1/3$. This also explains why the CB for SSS is highly likely to break down in
NS cores.

Figure 1: Sketch of the SSS near NS centers (panel (a)). The sign of
coefficient $D$ (see Eq. (30)) characterizes the possibility of
$s_{\textrm{c}}^{2}<s^{2}(\widehat{\varepsilon})$ (indication of a
crossover/soft core, panel (b)) or
$s^{2}_{\textrm{c}}>s^{2}(\widehat{\varepsilon})$ (the center EOS is stiffer
than its surroundings, panel (c)).

In this paper, we report results of our investigations into the remaining
issues listed above ((b)$\sim$(f)) in essentially the same framework as in our
earlier work [115]. When possible, we also compare our results with existing
ones by others in the literature. Our main contributions/findings include:

  1. (1)

The NS radius $R$ for a given dense matter EOS is shown to scale with
$\nu_{\textrm{c}}\equiv\varepsilon_{\textrm{c}}^{-1/2}[\widehat{P}_{\textrm{c}}/(1+3\widehat{P}_{\textrm{c}}^{2}+4\widehat{P}_{\textrm{c}})]^{1/2}$
along its M-R curve, extending the results of Ref. [115]. This scaling could
be used to reduce uncertainties of the core EOS for a certain NS considering
its mass/radius measurements. See FIG. 2 for illustrations of the R scaling.

  2. (2)

The SSS is shown to increase with $M_{\textrm{NS}}^{\max}$ under the
assumption that NSs with masses about
$M_{\textrm{NS}}^{\max}/M_{\odot}=1.3\mbox{$\sim$}2.3$ ($M_{\odot}=$ solar
mass) have similar radii about 12 km. As a result, the CB for SSS [116, 117]
is found to break down if $M_{\textrm{NS}}^{\max}\gtrsim 1.9M_{\odot}$. See
FIG. 5.

  3. (3)

Using the upper bound for the reduced central pressure
$\widehat{P}_{\textrm{c}}\lesssim 0.374$ to satisfy the causality requirement
and the universal correlations given in Ref. [115], a new causality boundary
for the NS M-R curve as $R_{\max}/\textrm{km}\gtrsim
4.73M_{\textrm{NS}}^{\max}/M_{\odot}+1.14$ is obtained. It is shown to be
consistent with observational data of several NSs albeit in certain tension
with predictions using a few empirical EOSs for NS matter in the literature.
As a direct corollary, the radius of the PSR J0952-0607 [118] with a mass
about $2.35M_{\odot}$ is found to be $\gtrsim$12.25 km. See FIG. 7 and FIG. 8.

  4. (4)

The ultimate energy density and pressure allowed in NSs are estimated, e.g.,
the existence of a $2.08M_{\odot}$ NS leads to
$\varepsilon_{\textrm{ult}}\lesssim 1.32\,\textrm{GeV}/\textrm{fm}^{3}$ and
$P_{\textrm{ult}}\lesssim 494\,\textrm{MeV}/\textrm{fm}^{3}$, which are found
to be consistent with several recent microscopic model analyses/calculations.
See FIG. 9.

  5. (5)

By adopting the empirical criteria
$\Theta=[(1/3-P/\varepsilon)^{2}+(s^{2}-P/\varepsilon)^{2}]^{1/2}\lesssim 0.2$
and $\gamma=\mathrm{d}\ln P/\mathrm{d}\ln\varepsilon\lesssim 1.75$ for
identifying conformality, we find it is unlikely to be realized in cores of
massive NSs. See TAB. 2, FIG. 16 and Eqs. (55)$\sim$(59).

  6. (6)

The $\widehat{P}_{\textrm{c}}$ is shown to be a key quantity determining the
onset of crossover at high densities (near the center). Specifically, the
probability of $s_{\textrm{c}}^{2}<s^{2}(\widehat{\varepsilon})$ for
$\widehat{\varepsilon}\lesssim 1$ (near the center) for the maximum-mass
configurations is found to be larger than about 50% for
$\widehat{P}_{\textrm{c}}\lesssim 0.3$. See FIG. 11, FIG. 12, FIG. 13 and FIG.
14.

  7. (7)

Continued with the last point, the existence of a peak in the derivative part
of $s^{2}$ as well as the non-existence of such peak in the non-derivative
part of $s^{2}$ are generally demonstrated, as introduced and discussed in
Ref. [119]. See Eqs. (50) and (52).

  8. (8)

By extending the upper bound $\widehat{P}_{\textrm{c}}\lesssim 0.374$ holding
at NS centers for the maximum-mass configurations, we demonstrated that the
reduced pressure $P/\varepsilon=\widehat{P}/\widehat{\varepsilon}$ is also
bounded from above to about 0.374 at finite distances near centers of stable
NSs. Such bound for $P/\varepsilon$ is a unique consequence of strong-field
gravity since it is much smaller in conventional (low-density) nuclear-physics
problems. As a direct corollary, the trace anomaly (TA) $\Delta\equiv
1/3-P/\varepsilon$ is bounded from below to about $\Delta\gtrsim-0.041$ near
the centers of massive NSs. See FIG. 17, Eqs. (47) and (49).

Here, our findings (1)$\sim$(5) concern properties at NS centers, while the
points (6)$\sim$(8) are extended to positions away from the centers (see FIG.
10 for a notation guide). A few empirical relations in the literature are also
explained using our method (e.g., FIG. 3 and FIG. 4). Since our results are
obtained without using any specific input EOS for NS matter, they provide
generally EOS-model independent insights on several interesting properties of
NSs especially near their cores where direct observational data is lacking.

The rest of this paper is organized as follows: in Section II we review
briefly the methods of Ref. [115] on directly extracting the NS central EOS
from observational data and then in Section III we estimate the central SS as
a function of $M_{\textrm{NS}}^{\max}$ using Eq. (1); we study the causality
boundary for the NS M-R curve in Section IV and the closely related
compactness parameter in Section V; and Section VI gives the estimate on the
ultimate (maximum) energy density as well as the pressure allowed in NS cores.
Section VII is devoted to the analysis on the possible occurrence of
continuous crossover in NS cores and in Section VIII and Section IX we
generalize the limit $\widehat{P}_{\textrm{c}}\lesssim 0.374$ which holds at
NS centers for the maximum-mass configurations to stable NSs along the M-R
curve either at centers or finite distances away. We then summarize this work
in Section X. In the appendices, we give relevant mathematical details of
expressions in the main text and estimates on certain quantities involved.

##  II A Brief Review: Neutron Star Mass and Radius Scalings from Perturbative
Solutions of the Scaled TOV Equations

For completeness and ease of our following discussions, we recall here the
main points of solving perturbatively the scaled TOV equations and the
resulting mass and radius scalings for NSs at the maximum-mass configuration.
More details of our approach can be found in Ref. [115]. We also discuss here
intuitively some general features of NSs from the scaled TOV equations that
are independent of the EOS-model and present the general radius scaling on the
NS mass-radius curve.

The TOV equations [120, 121, 122] describe the evolution of the pressure $P$
and mass $M$ as functions of the distance $r$ from the NS centers (adopting
$c=1$),

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}r}P=$
$\displaystyle-\frac{G(\varepsilon+P)(M+4\pi
r^{3}P)}{r^{2}(1-2GM/r)},~{}~{}\frac{\mathrm{d}}{\mathrm{d}r}M=4\pi
r^{2}\varepsilon.$ (2)

They are conventionally integrated from the center to surface (given a EOS for
NS matter) to give the NS mass $M_{\textrm{NS}}\equiv M$ and radius $R$
defined as the vanishing point of pressure (i.e., $P(R)=0$). We can rewrite
the TOV equations by scaling the NS mass by $W\equiv G^{-1}(4\pi
G\varepsilon_{\textrm{c}})^{-1/2}$ and the radius by $Q\equiv(4\pi
G\varepsilon_{\textrm{c}})^{-1/2}$ [115],

$\frac{\mathrm{d}}{\mathrm{d}\widehat{r}}\widehat{P}=-\frac{(\widehat{P}+\widehat{\varepsilon})(\widehat{r}^{3}\widehat{P}+\widehat{M})}{\widehat{r}^{2}-2\widehat{M}\widehat{r}},~{}~{}\frac{\mathrm{d}}{\mathrm{d}\widehat{r}}\widehat{M}=\widehat{r}^{2}\widehat{\varepsilon},$
(3)

here $\widehat{P}\equiv P/\varepsilon_{\textrm{c}}$ and
$\varepsilon_{\textrm{c}}$ is the central energy density. The above
dimensionless TOV equations (3) can be solved by expanding the (reduced)
energy density, pressure and mass as
$\widehat{\varepsilon}=1+\sum_{k=1}^{K}a_{k}\widehat{r}^{k}$,
$\widehat{P}=\widehat{P}_{\textrm{c}}+\sum_{k=1}^{K}b_{k}\widehat{r}^{k}$ and
$\widehat{M}=\sum_{k=1}^{K}c_{k}\widehat{r}^{k}$ with $\widehat{r}=r/Q$, here
$K$ is the effective truncation order of the polynomial expansions. As shown
in detail in Ref. [115], $c_{1}=c_{2}=0$ and $c_{k}=a_{k-3}/k$ for $k\geq 3$,
and the leading nonzero coefficient is
$b_{2}=-6^{-1}(1+3\widehat{P}_{\textrm{c}}^{2}+4\widehat{P}_{\textrm{c}})$.
Truncating the pressure to order $b_{2}$, namely
$\widehat{P}_{\textrm{c}}+b_{2}\widehat{R}^{2}=0$, gives the reduced radius
$\widehat{R}\sim[\widehat{P}_{\textrm{c}}/(1+3\widehat{P}_{\textrm{c}}^{2}+4\widehat{P}_{\textrm{c}})]^{1/2}$
and the physical radius $R\sim\widehat{R}Q$ is found to scale according to
[115]

$R\sim\nu_{\textrm{c}}\equiv\frac{\widehat{P}_{\textrm{c}}^{1/2}}{\sqrt{\varepsilon_{\textrm{c}}}}\left(\frac{1}{1+3\widehat{P}_{\textrm{c}}^{2}+4\widehat{P}_{\textrm{c}}}\right)^{1/2}.$
(4)

The NS mass
$M_{\textrm{NS}}\sim\widehat{R}^{3}/\sqrt{\varepsilon_{\textrm{c}}}$ is found
to scale as [115]

$M_{\textrm{NS}}\sim\Gamma_{\textrm{c}}\equiv\frac{\widehat{P}_{\textrm{c}}^{3/2}}{\sqrt{\varepsilon_{\textrm{c}}}}\left(\frac{1}{1+3\widehat{P}_{\textrm{c}}^{2}+4\widehat{P}_{\textrm{c}}}\right)^{3/2}.$
(5)

In Ref. [115], the correlations $M_{\textrm{NS}}^{\max}$-$\Gamma_{\textrm{c}}$
and $R_{\max}$-$\nu_{\textrm{c}}$ are then verified using 87 widely used
phenomenological and 17 microscopic EOSs, here $R_{\max}$ is the corresponding
radius at $M_{\textrm{NS}}^{\max}$. Specifically [115],

$\displaystyle R_{\max}/\textrm{km}$ $\displaystyle\approx 1.05\times
10^{3}\nu_{\textrm{c}}+0.64,$ (6) $\displaystyle
M_{\textrm{NS}}^{\max}/M_{\odot}$ $\displaystyle\approx 0.173\times
10^{4}\Gamma_{\textrm{c}}-0.106.$ (7)

The SSS of Eq. (1) is obtained via
$\mathrm{d}M_{\textrm{NS}}/\mathrm{d}\varepsilon_{\textrm{c}}=0$ at the
maximum-mass configuration on the M-R curve.

Figure 2: Correlation between the NS radii $R$ and the factor
$\nu_{\textrm{c}}$ obtained from the non-relativistic EDFs (of the Gogny-like
and Skyrme types) and the nonlinear RMF model EOSs. The solid symbol is for
the maximum-mass configuration of each model (corresponding to the solid
symbol shown in the inset).

The scalings of (4) and (5) could be understood intuitively by combining the
quantum-degeneracy of nuclear pressure and the self-gravitating nature of NSs.
Since NSs are self-gravitating systems [123, 124, 125, 126, 127, 128, 129,
130, 131, 132], one expects that a larger energy density $\varepsilon$
corresponds to a smaller radius $R$. By temporally neglecting the general
relativistic effects [133], one has
$\mathrm{d}P/\mathrm{d}r=-GM\varepsilon/r^{2}$ and thus $P\sim-
GM\varepsilon/r\sim Gr^{2}\varepsilon^{2}$ where $\mathrm{d}M/\mathrm{d}r=4\pi
r^{2}\varepsilon$ or $M\sim r^{3}\varepsilon$ is used. Thus $r\sim
1/\sqrt{G\varepsilon}$ because $P$ and $\varepsilon$ have the same dimension.
On the other hand, NSs are supported mainly by the neutron degenerate
pressure, a larger pressure $P$ naturally leads to a larger radius $R$, i.e.,
$R\sim P^{\phi}$ where $\phi>0$. In order to infer the value of $\phi$, one
notices that $P\sim\int M\varepsilon/r^{2}\mathrm{d}r$ as a function of $r$ is
even, since the $\varepsilon$ is even of $r$ and thus the mass $M$ is odd of
$r$ because $M\sim\int r^{2}\varepsilon\mathrm{d}r$. Consequently,
$P\sim\textrm{const}.+Br^{2}+\cdots$, from which one obtains $\phi=1/2$ and
$R\sim P^{1/2}+\cdots$. The absence of linear term ($\propto r$) in the
expansion of $P$ over $r$ could also be understood by the boundary condition
$\mathrm{d}P/\mathrm{d}r=0$ at $r=0$ (pressure cannot have a cusp-like
singularity). By combining the self-gravitating and quantum-degenerate nature
of NSs, we have $R\sim(P/\varepsilon)^{1/2}/\sqrt{G\varepsilon}+\cdots$, since
$P/\varepsilon$ is the relevant dimensionless quantity from combining $P$ and
$\varepsilon$. The NS M-R relation is obtained by integrating the TOV
equations from a given central energy density $\varepsilon_{\textrm{c}}$, thus
$R\sim\widehat{P}_{\textrm{c}}^{1/2}/\sqrt{G\varepsilon_{\textrm{c}}}+\cdots\sim\widehat{P}_{\textrm{c}}^{1/2}/\sqrt{\varepsilon_{\textrm{c}}}+\cdots$
where $\widehat{P}_{\textrm{c}}=P_{\textrm{c}}/\varepsilon_{\textrm{c}}$.
Considering general relativistic effects, the scaling for $R$ will be modified
to
$R\sim\widehat{P}_{\textrm{c}}^{1/2}/\sqrt{\varepsilon_{\textrm{c}}}\cdot\vartheta(\widehat{P}_{\textrm{c}})$
with the correction $\vartheta(\widehat{P}_{\textrm{c}})$ to be revealed by
analyzing structures of the TOV equations [115]. In particular, Eq. (4) tells
that
$\vartheta(\widehat{P}_{\textrm{c}})=(1+3\widehat{P}_{\textrm{c}}^{2}+4\widehat{P}_{\textrm{c}})^{-1/2}<1$,
i.e., the strong gravity in GR than Newtonian’s reduces the NS radius. Similar
arguments give for the mass as $M_{\textrm{NS}}\sim
R^{3}\varepsilon_{\textrm{c}}\sim\widehat{P}^{3/2}/G^{3/2}\sqrt{\varepsilon_{\textrm{c}}}+\cdots\sim\widehat{P}_{\textrm{c}}^{3/2}/\sqrt{\varepsilon_{\textrm{c}}}+\cdots$.

At this point, it is necessary to point out the potential degrees of
uncertainty inherent in our results, due to the pertubative nature of our
approach. While some of them are general and robust, e.g., the correlations
between $M_{\textrm{NS}}^{\max}$ and $\Gamma_{\textrm{c}}$ and that between
$R_{\max}$ and $\nu_{\textrm{c}}$ [115] and the consequent prediction on the
causality boundary of FIG. 7, the compactness constraint of FIG. 8 and the
ultimate energy density and pressure allowed in NSs of FIG. 9, others may
contain uncertainties. For example, when estimating the possible peak of
$s^{2}$ as a function of $\widehat{r}$ (see Eq. (39)), we expand the
$s^{2}(\widehat{r})$ to order $\widehat{r}^{4}$ which neglects the higher-
order contributions from $\widehat{r}^{6}$ and $\widehat{r}^{8}$, etc., and
these higher-order terms may induce quantitative effects on the estimate for
the peak. However, when investigating whether there exists a peak in
$s^{2}(\widehat{r})$, the estimate on the sign of the coefficient of
$\widehat{r}^{2}$ in the expansion of $s^{2}(\widehat{r})$ is enough, and in
this sense the result avoids suffering from contributions from the higher-
order terms. Similarly, when bounding the ratio $P/\varepsilon$ in Section
VIII, the conclusion is effective only near the NS cores and whether it is
still effective at places far from the center needs further analysis by
including higher-order terms in $\widehat{r}$. When necessary, we may point
out the relevant uncertainties and/or the limitations of our results, see,
e.g., the last second paragraph of Section VII.

Besides the correlation between $M_{\textrm{NS}}^{\max}$ and
$\Gamma_{\textrm{c}}$ discussed in Ref. [115], the NS radius for a given dense
matter EOS is also found to be strongly correlated with $\nu_{\textrm{c}}$
along the M-R curve. We show in FIG. 2 the dependence of radius $R$ on
$\nu_{\textrm{c}}$ using three kinds of EOSs: the covariant relativistic mean-
field (RMF) models [134] and the non-relativistic energy density functionals
(EDFs) of the Gogny-like type [135] as well as the conventional Skyrme-type
[136, 138, 137]. Very good correlations between $R$ and $\nu_{\textrm{c}}$ are
found in all cases as indicated by the very similar slopes (444, 454 and 462
in unit of $[R]/[10^{3}\nu_{\textrm{c}}]$ shown by the captions near the
lines), while the intercept has certain model-dependence (considering the two
classes). A cutoff of $M_{\textrm{NS}}\gtrsim 1.2M_{\odot}$ is necessary to
mitigate potential influences of uncertainties in modeling the crust EOS [139,
140, 141] for low-mass NSs. The radii $R_{\max}$ are also shown using solid
symbols ($M_{\textrm{NS}}^{\max}/M_{\odot}\approx 2.05$ in the Skyrme model
using green diamond). Points on the left side of $R_{\max}$ correspond to
unstable NS configurations (i.e., for such NSs one has
$\mathrm{d}M_{\textrm{NS}}/\mathrm{d}\varepsilon_{\textrm{c}}<0$). The
$R$-$\nu_{\textrm{c}}$ correlation could be used to reduce uncertainties of
the central EOS for a NS with certain mass via its radius
constraints/measurements within a specific EOS model.

Figure 3: Correlation between $\rho_{\textrm{c}}/\rho_{\textrm{sat}}$ and
radius $R_{\max}$ adopting two reference values for $\widehat{P}_{\textrm{c}}$
(0.16 and 0.24), here $\rho_{\textrm{sat}}\approx 0.16\,\textrm{fm}^{-3}$ is
used. An empirical prediction on
$\rho_{\textrm{c}}/\rho_{\textrm{sat}}$-$R_{\max}$ from Ref. [143] (orange
solid contour) is shown for a comparison.

By using the scaling in Eq. (6) we can derive a relation between NS central
baryon density $\rho_{\textrm{c}}$ and the radius $R_{\max}$. The result is,

$\frac{\rho_{\textrm{c}}}{\rho_{\textrm{sat}}}\approx\frac{7.35\times
10^{3}\widehat{P}_{\textrm{c}}}{1+3\widehat{P}_{\textrm{c}}^{2}+4\widehat{P}_{\textrm{c}}}\left(\frac{R_{\max}}{\textrm{km}}-0.64\right)^{-2},$
(8)

where $\rho_{\textrm{sat}}\approx 0.16\,\textrm{fm}^{-3}$ [142] is the nuclear
saturation density and the approximation
$\rho_{\textrm{c}}\approx\varepsilon_{\textrm{c}}/M_{\textrm{N}}$ is adopted
with $M_{\textrm{N}}\approx 939\,\textrm{MeV}$ the static nucleon mass. Shown
in FIG. 3 is the relation between $\rho_{\textrm{c}}/\rho_{\textrm{sat}}$ and
the radius $R_{\max}$ using two values of $\widehat{P}_{\textrm{c}}$ (0.16 and
0.24). For a comparison, an algorithmic/empirical prediction from Ref. [143]
is also shown (orange solid contour). It is seen that the overall correlation
between $\rho_{\textrm{c}}/\rho_{\textrm{sat}}$ and $R_{\max}$ from our
analyses here is consistent with that from Ref. [143]. Moreover, one can see
from Eq. (8) that roughly $\rho_{\textrm{c}}/\rho_{\textrm{sat}}\sim
R_{\max}^{-2}\cdot[1+\mbox{corrections of }R_{\max}^{-1}]$ once a
$\widehat{P}_{\textrm{c}}$ is specified. The correlation shown in FIG. 3 may
help us to estimate the maximum central (baryon) density when future radius
measurements of very massive NSs are available [144]. Moreover, the above
correlation is obtained under the approximation
$\varepsilon_{\textrm{c}}\approx M_{\textrm{N}}\rho_{\textrm{c}}$ and
principally one needs to use $E(\rho_{\textrm{c}})+M_{\textrm{N}}$ to replace
$M_{\textrm{N}}$ where $E(\rho_{\textrm{c}})$ defined as the energy per
nucleon at $\rho_{\textrm{c}}$ has its own uncertainties. Doing so will shift
down the lines in FIG. 3 as $E(\rho_{\textrm{c}})>0$.

Figure 4: A comparison of the radius scalings ($R$ vs reduced pressure
$\widehat{P}$): the empirical power law $R\sim P^{\delta}$ with $\delta\approx
0.23\mbox{$\sim$}0.26$ from Ref. [13], the Newtonian scaling
$\widehat{R}\sim\widehat{P}^{1/2}$ and our full scaling
$\widehat{R}\sim[\widehat{P}/(1+3\widehat{P}^{2}+4\widehat{P})]^{1/2}$.

It is appropriate at this point to compare our radius-pressure scaling with
existing ones in the literature in FIG. 4. The grey band shows the empirical
power law that the NS radii depend on the pressure at 1-2 times
$\rho_{\textrm{sat}}$, i.e., $R\sim P^{\delta}$ with $\delta\approx
0.23\mbox{$\sim$}0.26$ [13]. One finds that the Newtonian prediction
$\widehat{R}\sim\widehat{P}^{1/2}$ (by neglecting the general-relativistic
correction $3\widehat{P}^{2}+4\widehat{P}$) indicated by the blue dashed line
deviates significantly from the empirical power law $\widehat{P}^{\delta}$. On
the other hand, our full scaling
$\widehat{R}\sim[\widehat{P}/(1+3\widehat{P}^{2}+4\widehat{P})]^{1/2}$
(magenta line) is rather consistent with the empirical scaling, especially for
$0.1\lesssim\widehat{P}\lesssim 0.3$ (which is a reasonable region for
$P/\varepsilon$ in NS cores).

##  III Speed of Sound and Violation of the Conformal Bound in Massive Neutron
Stars

The dependence of $s_{\textrm{c}}^{2}$ on $\widehat{P}_{\textrm{c}}$ could be
straightforwardly transformed into its dependence on
$M_{\textrm{NS}}^{\max}/M_{\odot}$. In the following, we assume that NSs with
masses $\textrm{X}_{-0.07}^{+0.07}M_{\odot}$ where
$\textrm{X}=1.3\mbox{$\sim$}2.3$ have similar radii as indicated by the NICER
observations to extract the $\nu_{\textrm{c}}$ from Eq. (6). Specifically,
NICER found that the radius of PSR J0740+6620 (mass $\approx
2.08_{-0.07}^{+0.07}M_{\odot}$) is about $12.39_{-0.98}^{+1.30}\,\textrm{km}$
[70] and that of the canonical PSR J0030+0451 (mass $\approx
1.34_{-0.16}^{+0.15}M_{\odot}$) is about $12.71_{-1.19}^{+1.14}\,\textrm{km}$
[73] (see also Ref. [74]). They are very similar, indicating that the NS
radius is less sensitive to the central EOS [115]. Therefore, we obtain
$10^{3}\nu_{\textrm{c}}\approx 11.2\,\textrm{fm}^{3/2}/\textrm{MeV}^{1/2}$ by
using $R_{\max}\approx 12.39_{-0.98}^{+1.30}\,\textrm{km}$ [70].
Perturbatively, we have for the SSS at NS center from expanding Eq. (1) as,

$\displaystyle s_{\textrm{c}}^{2}$
$\displaystyle\approx\frac{4}{3}\widehat{P}_{\textrm{c}}\left(1+\widehat{P}_{\textrm{c}}+\frac{3}{2}\widehat{P}_{\textrm{c}}^{2}+3\widehat{P}_{\textrm{c}}^{3}\right)+\mathcal{O}\left(\widehat{P}_{\textrm{c}}^{5}\right)$
(9)
$\displaystyle\approx\frac{4}{3}\beta\left(1+5\beta+\frac{57}{2}\beta^{2}+175\beta^{3}\right)+\mathcal{O}\left(\beta^{5}\right),$
(10)

where $\beta\approx 0.052(M_{\textrm{NS}}^{\max}/M_{\odot}+0.106)\ll 1$.
Keeping only the leading-order term of (9) gives the Newtonian prediction
$s_{\textrm{c}}^{2}=4\widehat{P}_{\textrm{c}}/3$ from which one infers
$\widehat{P}_{\textrm{c}}\lesssim 3/4$. The general-relativistic contributions
in the TOV equations have an effect about 100% on the upper limit for
$\widehat{P}_{\textrm{c}}$. Essentially, the $s_{\textrm{c}}^{2}$ increases as
$M_{\textrm{NS}}^{\max}$ increases as shown in FIG. 5. Results of FIG. 5
provide us a straightforward way to infer the central sound speed
$s_{\textrm{c}}^{2}$ once the radii/masses are known (measured). It also
implies that more compact NSs have larger $s_{\textrm{c}}^{2}$, e.g., NSs of
masses $1.5M_{\odot}$ and $1.4M_{\odot}$ with radii 9.9$\sim$11.2 km [145] and
$11.0_{-0.6}^{+0.9}\,\textrm{km}$ [146] have $s_{\textrm{c}}^{2}\approx
0.31_{-0.07}^{+0.07}$ and $s_{\textrm{c}}^{2}\approx 0.23_{-0.05}^{+0.03}$,
respectively, shown as red hexagons in FIG. 5. While for PSR J0740+6620 as
indicated with the green diamond, we have $s_{\textrm{c}}^{2}\approx
0.45_{-0.18}^{+0.14}$ using $\widehat{P}_{\textrm{c}}\approx
0.24_{-0.07}^{+0.05}$ [115].

Figure 5: Dependence of $s_{\textrm{c}}^{2}$ on the maximum mass
$M_{\textrm{NS}}^{\max}/M_{\odot}$ (of a given EOS for NS matter), the CB on
SSS is indicated by the dashed orange line. The inset plots the
$s_{\textrm{c}}^{2}$ as a function of $\widehat{P}_{\textrm{c}}$, where the
lightblue curve is allowed by conditions of stability and causality ($0\leq
s_{\textrm{c}}^{2}\leq 1$).

In FIG. 5, the CB for SSS [116, 117], namely $s^{2}=1/3$ is also indicated (by
the dashed orange line). Our result on $s_{\textrm{c}}^{2}$ implies that the
CB for SSS tends to break down for $M_{\textrm{NS}}^{\max}/M_{\odot}\gtrsim
1.89$. The latter is very close to the predicted most probable critical
maximum-mass about $M_{\textrm{NS}}^{\max}/M_{\odot}\approx 1.87$ given in
Ref. [147], above which the CB is likely to break down (see their FIG. 1 and
FIG. 2). See also Ref. [154], which predicted that the largest mass of NSs
being consistent with the CB should be $\lesssim 1.99M_{\odot}$. Equivalently,
the CB $s_{\textrm{c}}^{2}\leq 1/3$ is broken for
$\widehat{P}_{\textrm{c}}\gtrsim 0.195$ from Eq. (1). The finding here is
consistent with previous studies on the same issue [153, 149, 150, 154, 151,
152, 155, 156, 157, 148], indicating that Eq. (1) grasps the main features of
the SS. However, we emphasize that our method adopts no specific model for NS
matter EOS. We also make no assumptions about the composition of NSs
(nucleons, hyperons and/or quarks), and therefore provides a EOS-model
independent way for investigating the SS in NSs. Since the NS EOSs are often
softened considering the exotic components such as hyperons, our formula (1)
indicates that the CB in these NSs containing non-nucleonic particles is
relatively easier to be obeyed, see, e.g., Refs. [158, 159, 160] on related
issues.

To avoid causing confusion, we emphasize here that although NSs with lower
maximum-mass have $s_{\textrm{c}}^{2}\lesssim 1/3$, it does not imply that
these NSs exhibit conformal symmetry (in their cores) since the latter emerges
at extremely high energies where the strong force coupling constant
diminishes, leading to the “asymptotic freedom” of quarks and gluons. This
symmetry is broken as the energy scale lowers or the theory is probed at
larger scales. At the energy scales relevant to NS cores, the presence of
massive particles and strong gravitational fields disrupt conformal symmetry.
The latter being not exact even in quark-gluon plasma is generally expected to
break down as the system transitions to a hadronic phase. Thus, the
observation of $s_{\textrm{c}}^{2}\lesssim 1/3$ in low-mass NSs should not be
interpreted to imply that the conformal symmetry is reached there [161].

##  IV Causality Boundary for Neutron Star Mass-radius Curve

The correlations of $M_{\textrm{NS}}^{\max}$-$\Gamma_{\textrm{c}}$ and
$R_{\max}$-$\nu_{\textrm{c}}$ also enable us to localize the causality
boundary of the NS M-R curve. Combining Eqs. (6) and (7) leads to,

$\frac{M_{\textrm{NS}}^{\max}}{M_{\odot}}\approx\frac{1.65\widehat{P}_{\textrm{c}}}{1+3\widehat{P}_{\textrm{c}}^{2}+4\widehat{P}_{\textrm{c}}}\left(\frac{R_{\max}}{\textrm{km}}-0.64\right)-0.106.$
(11)

The in-front coefficient $f\approx
1.65\widehat{P}_{\textrm{c}}/(1+3\widehat{P}_{\textrm{c}}^{2}+4\widehat{P}_{\textrm{c}})$
is a slow-varying function of $\widehat{P}_{\textrm{c}}$, and it essentially
explains the conventional quasi-linear correlation between
$M_{\textrm{NS}}^{\max}/M_{\odot}$ and $R_{\max}$ from model calculations. In
FIG. 6, we show the scatters of $M_{\textrm{NS}}^{\max}$ and $R_{\max}$ using
the RMF approaches and the non-relativistic EDFs of Gogny-like and Skyrme
types together with a few microscopic EOSs [115]. The fitting lines of Eq.
(11) adopt four fiducial values for $\widehat{P}_{\textrm{c}}$ (captioned near
the lines with the same color), from which one finds that for
$\widehat{P}_{\textrm{c}}\gtrsim 0.3$ the expression (11) could reasonably
describe the EOS samples. However, obvious dispersions can be seen for these
scatters. In addition, certain EOSs are unfavored by the causal limit
$\widehat{P}_{\textrm{c}}\lesssim 0.374$ (solid black line, corresponding to
$s_{\textrm{c}}^{2}\leq 1$ [115]), see FIG. 7 for further discussions.

Figure 6: Scatters of $M_{\textrm{NS}}^{\max}$ and $R_{\max}$ with the
covariant RMF approaches, the non-relativistic EDFs of Gogny-like and Skyrme
types and several extra microscopic EOSs [115]. The captions near each lines
are the fiducial values for $\widehat{P}_{\textrm{c}}$.  Figure 7: Causality
boundaries for NS M-R curve, here the black solid line is based on inequality
(12) (see the corresponding black line in FIG. 6). The prediction using the
trace anomaly $\Delta\equiv 1/3-P/\varepsilon\geq 0$ (blue dashed line) as
well as the previous constraint $R_{\max}/\textrm{km}\gtrsim
4.51M_{\textrm{NS}}^{\max}/M_{\odot}$ (grey dotted line) are also shown for
comparison. The M-R constraints for several NSs (magenta/green/lightblue/red
contours), M-R curves from a few typical empirical dense matter EOSs (grey
solid lines), and the constraints for the M-R relation based on events
GW170817 (chocolate/pink bands) and GW190425 (cyan solid contour) are also
shown. See text for details.

Since the factor
$(1+3\widehat{P}_{\textrm{c}}^{2}+4\widehat{P}_{\textrm{c}})/\widehat{P}_{\textrm{c}}$
takes its minimum about 7.80 at $\widehat{P}_{\textrm{c}}\approx 0.374$, we
have from Eq. (11) that,

${R_{\max}}/{\textrm{km}}\gtrsim
4.73{M_{\textrm{NS}}^{\max}}/{M_{\odot}}+1.14.$ (12)

The inequality (12) is obtained for the maximum-mass configuration via the
causal condition $s_{\textrm{c}}^{2}\leq 1$, it therefore gives a causality
boundary for the NS M-R curve. In FIG. 7, we show the causality boundary of
(12) (black solid line) together with the M-R constraints for a few NSs
including the J0740+6620 and J0030+0451 [69, 70, 71, 72, 73, 74]
(red/lightblue contours), the 4U 1820-30 and EXO 1745-248 [162] (green/magenta
contours), the recently reported (“black widow”) PSR J0952-0607 with a mass
about 2.35$M_{\odot}$ [118] (pink hatched band), GW190814’s secondary
component with a mass of $2.59^{+0.08}_{-0.09}M_{\odot}$ [163] (cyan hatched
band) and the M-R curves from several empirical dense matter EOSs [164, 167,
166, 8, 27, 168, 171, 170, 165, 172, 173, 174, 169, 175] (grey solid lines).
The violet and the chocolate bands are obtained from the GW170817 event using
the universal relations (URs) and the spectral EOS (specEOS) approaches [47],
respectively, while the cyan solid contour is based on the GW190425 event for
a compact binary coalescence with total mass about 3.4M⊙ [176]. Additionally,
the grey dotted line (marked as “LPMY”) was empirically given by Ref. [177] as
$R_{\max}/\textrm{km}\gtrsim 4.51M_{\textrm{NS}}^{\max}/M_{\odot}$, by
considering nuclear EOS effects. Recently, the causality boundary for NS M-R
curve was localized using an analysis of TA [119] via the condition
$\Delta\geq 0$ (blue dashed line), where $\Delta\equiv 1/3-P/\varepsilon$ was
suggested to measure the TA (see Section VIII for more discussions on $\Delta$
and related issues).

At this point, it is necessary to discuss in some more details the fundamental
difference between the apparent condition $P/\varepsilon\leq 1$ from the
principle of causality and $s^{2}\leq 1$, with the latter being quite relevant
in NSs. The physical origin of this distinction can be traced back to the
nontrivial/nonlinear dependence of the $s_{\textrm{c}}^{2}$ on
$\widehat{P}_{\textrm{c}}$, as shown clearly in Eq. (1) which was revealed by
the intrinsic structures of the TOV equations [115] (similar clarifications
were given recently in Ref. [119] using the TA). Actually, even an “ultra-
stiff” EOS $P=\varepsilon$ may introduce probable tension/conflict with the
observational data of NSs, e.g., see and compare the upper grey dashed line
(marked by “$P=\varepsilon$”) and the red solid contour for PSR J0740+6620 (at
a 68% confidence level). Specifically, one can solve analytically the TOV
equations adopting the linear EOS $P=w\varepsilon$ (for which
$s^{2}=P/\varepsilon=w$ and therefore $P/\varepsilon\leq 1$ and $s^{2}\leq 1$
are totally equivalent) to give the M-R relation as (adopting $c=G=1$),

$M_{\textrm{NS}}^{(w)}=\frac{2w}{1+6w+w^{2}}R.$ (13)

Consequently, $M_{\textrm{NS}}^{(w=1)}=R/4$ is obtained, which is smaller than
both the Schwarzschild relation $M_{\textrm{Schw}}=R/2$ and the Buchdahl’s
prediction $M_{\textrm{Buch}}=4R/9$ [178], see the tan and deep-tan bands for
the latter two relations in FIG. 7. On the other hand, compared with the
observational data, the ultra-relativistic Fermi gas EOS (with
$P=\varepsilon/3$) leads to quite an inconsistent causality boundary (shown as
the lower grey dashed line marked by “$P=\varepsilon/3$”), i.e.,
$M_{\textrm{NS}}^{(w=1/3)}=3R/14$ [179]. The central SSS is $7/9$ if
$\widehat{P}_{\textrm{c}}=1/3$ is inserted into Eq. (1).

Interestingly, our result (12) puts a more stringent constraint on the M-R
curve than those through the TA analysis [119] and the one from Ref. [177].
Moreover, it is very consistent with the observational data while excluding a
few celebrated empirical dense matter EOSs as indicated in FIG. 7. In
particular, the observational boundary of M-R data for PSR J0740+6620 [69, 70,
71, 72] and those from the events GW170817 [47] and GW190425 [176] are
basically consistent with the relation (12). Therefore, it is important to
check the new causality boundary (12) as more observational data on NS
masses/radii or GW events become available (see similar discussions given in
Ref. [119] and their FIG.3). As an illustration of applying the new causal
boundary on the mass-radius curve, we give a lower limit for the radius as
$R_{\max}\gtrsim 12.25\,\textrm{km}$ for the recently identified PSR
J0952-0607 [118] from (12), while for the canonical NS the radius is bounded
from below to about 7.77 km (shown as the cyan solid diamond in FIG. 7).
Inversely, if the radii of NSs were known (constrained), e.g.,
$R_{\max}\approx 12.39\,\textrm{km}$ (for PSR J0740+6620 [70] from NICER),
using the causal limit $\widehat{P}_{\textrm{c}}\lesssim 0.374$ gives an upper
limit $M_{\textrm{NS}}^{\max}\lesssim 2.38M_{\odot}$. Similarly if
$R_{\max}\approx 11.41\,\textrm{km}$ or $13.69\,\textrm{km}$ (the lower/upper
constraints for radius of PSR J0740+6620 [70]) is adopted, one has
$M_{\textrm{NS}}^{\max}\lesssim 2.17M_{\odot}$ or
$M_{\textrm{NS}}^{\max}\lesssim 2.65M_{\odot}$, respectively.

##  V Constraints on Neutron Star Compactness

An important quantity connected with the causality boundary is the NS
compactness (parameter) defined as $\xi=M_{\textrm{NS}}/R$ (adopting $c=G=1$).
For the maximum-mass configuration on the NS M-R curve, one obtains the upper
bound for $\xi_{\max}={M_{\textrm{NS}}^{\max}}/{R_{\max}}$ by neglecting the
intercepts $0.106$ and $0.64$ in Eq. (7) and Eq. (6) respectively, as

$\xi_{\max}\leq\frac{0.173\times 10^{4}\Gamma_{\textrm{c}}}{1.05\times
10^{3}\nu_{\textrm{c}}}\left(\frac{M_{\odot}}{\textrm{km}}\right)\approx\frac{2.44\widehat{P}_{\textrm{c}}}{1+3\widehat{P}_{\textrm{c}}^{2}+4\widehat{P}_{\textrm{c}}}\equiv\xi_{\max}^{\textrm{(up)}}.$
(14)

Therefore, $\xi_{\max}\lesssim 0.313$ by taking
$\widehat{P}_{\textrm{c}}\approx 0.374$, which is about $30\%$ smaller than
the Buchdahl limit $4/9\approx 0.444$.

Figure 8: Compactness parameter $\xi_{\max}$ for the maximum-mass
configuration as a function of central SSS $s_{\textrm{c}}^{2}$. Several other
constraints/observational samples (see the text for details) are also shown
for comparisons.

In FIG. 8, we plot the compactness parameter $\xi_{\max}$ of Eq. (14) as a
function of the central SSS $s_{\textrm{c}}^{2}$ (which depends on the
$\widehat{P}_{\textrm{c}}$ in the form of Eq. (1)), where the compactness for
PSR J0740-6620 about $0.217\mbox{$\sim$}0.279$ (with the central value about
0.248) is also shown by the lightblue hatched band (and the purple shallow
circle) directly from its observational mass about
$2.08_{-0.07}^{+0.07}M_{\odot}$ [69] and radius about
$12.39_{-0.98}^{+1.30}\,\textrm{km}$ [70]. The compactness about
$0.163_{-0.026}^{+0.019}$ for PSR J0030+0451 [74] is shown by the light-green
diamond (the SSS is obtained from its mass and radius using Eqs. (7), (6) and
(1)). A similar result for NS EXO 1745-248 [145] is shown by the magenta
hexagon. We notice that neglecting the intercepts $0.106$ and $0.64$ in Eq.
(7) and Eq. (6) slightly overestimates the compactness. Actually, by rewriting
(12) we can obtain

$\xi_{\max}\lesssim 0.313\cdot\left(1-{1.14\,\textrm{km}}/{R_{\max}}\right),$
(15)

where $1.14/R_{\max}\ll 1$ could be treated as a correction to Eq. (14).
Considering a typical NS radius $R_{\max}\approx 12\,\textrm{km}$ (with a
reduction about 10% on $\xi_{\max}$), it then leads to $\xi_{\max}\lesssim
0.283$ (grey dash-dotted line). In this sense, the $\xi_{\max}\lesssim 0.313$
from Eq. (14) provides an upper limit for the NS compactness. An early
constraint on the compactness $\xi$ about $\xi\lesssim 0.353$ [127] is
indicated by the dashed red line (with a row of arrows). Furthermore, the
constraint on NS compactness can also be obtained from the gravitational
redshift $z$ defined by $1+z=(1-2\xi)^{-1/2}$. Interestingly, a new constraint
on $z$ was derived very recently from comparing X-ray burst model simulations
with observational data for GS 1826-24 using newly measured atomic masses
around the rp-process waiting-point nucleus ${}^{64}\textrm{Ge}$ [180]. The
resulting NS compactness is about $0.244\mbox{$\sim$}0.342$ with the most
probable value 0.27 at 95% confidence level (shown by the tan band). In
addition, the lower limit of the tan band and the dash-dotted grey line
together imply that $s_{\textrm{c}}^{2}\gtrsim 0.4$, implying the CB of 1/3
tends to break down. Compared to all of the above constrains available on NS
compactness, ours is probably the most stringent one so far.

We notice by passing that the compactness parameter is closely related to some
bulk properties of NSs such as the moment of inertial $I$ and/or the NS
binding energy $E_{\textrm{b}}$ [13, 14, 15, 141, 181, 183, 182, 184].
Therefore, the constraint on $\xi$ may naturally lead to some constraints on
these quantities. For example, via the empirical expression for the moment of
inertial as $I/M_{\textrm{NS}}R^{2}\approx 0.237(1+2.84\xi^{2}+18.9\xi^{4})$
given by Ref. [182], one may obtain that $I/M_{\textrm{NS}}R^{2}\lesssim 0.49$
using $\xi\lesssim 0.313$ or $I/M_{\textrm{NS}}R^{2}\lesssim 0.46$ using
$\xi\lesssim 0.283$. Similarly, Ref. [13] approximated the $E_{\textrm{b}}$ as
$E_{\textrm{b}}/M_{\textrm{NS}}\approx 0.6\xi/(1-\xi/2)$, and therefore
$E_{\textrm{b}}/M_{\textrm{NS}}\lesssim 0.22$ and
$E_{\textrm{b}}/M_{\textrm{NS}}\lesssim 0.20$ using $\xi\lesssim 0.313$ and
$\xi\lesssim 0.283$, respectively. We will not discuss these quantities
further in the current work since they could be treated as direct corollaries
of the constraints on $\xi$.

##  VI The Ultimate Energy Density and Pressure in Neutron Stars

Actually, the correlation $M_{\textrm{NS}}^{\max}$-$\Gamma_{\textrm{c}}$ alone
is already useful for inferring the ultimate (maximum) energy density
$\varepsilon_{\textrm{ult}}$ as well as the ultimate pressure
$P_{\textrm{ult}}$ allowed in NSs before they collapse into black holes. By
rewriting Eq. (7) and considering the casual limit
$\widehat{P}_{\textrm{c}}\lesssim 0.374$, we obtain

$\varepsilon_{\textrm{c}}\leq\varepsilon_{\textrm{ult}}\equiv
6.32\left({M_{\textrm{NS}}^{\max}}/{M_{\odot}}+0.106\right)^{-2}\,\textrm{GeV}/\textrm{fm}^{3}.$
(16)

Using the constraint $\widehat{P}_{\textrm{c}}\lesssim 0.374$ once again
gives,

$P_{\textrm{c}}\leq P_{\textrm{ult}}\equiv
2.36\left({M_{\textrm{NS}}^{\max}}/{M_{\odot}}+0.106\right)^{-2}\,\textrm{GeV}/\textrm{fm}^{3}.$
(17)

We note that the correlation between $R_{\max}$ and $\nu_{\textrm{c}}$ of Eq.
(6) is not used in establishing these limits.

Figure 9: The ultimate energy density $\varepsilon_{\textrm{ult}}$ and the pressure $P_{\textrm{ult}}$ (inset) allowed in NS cores as functions of the maximum-mass $M_{\textrm{NS}}^{\max}$ of NS M-R relation. The previous predictions on $\varepsilon_{\textrm{ult}}$ and $P_{\textrm{ult}}$ from Ref. [131] are also shown for comparison (black dashed line). See text for details.  $M_{\textrm{NS}}^{\max}/M_{\odot}$ | $\varepsilon_{\textrm{ult}}$ | $P_{\textrm{ult}}$ | Ref. | Eq. (16) | Eq. (17)
---|---|---|---|---|---
$\gtrsim 1.97$ | / | $570_{-320}^{+320}$ | [146] | 1.47 | 548
$\gtrsim 2.0$ | 1.46 | / | [154] | 1.42 | 533
$\gtrsim 2.01$ | $1.18_{-0.17}^{+0.17}$ | / | [187] | 1.41 | 528
$\gtrsim 2.2$ | $1.16_{-0.10}^{+0.11}$ | / | [188] | 1.19 | 444
TAB. 1: Summary of the constraints on $\varepsilon_{\textrm{ult}}$ and
$P_{\textrm{ult}}$ existing in the literature and from Eqs. (16) and (17),
here $\varepsilon_{\textrm{ult}}$ and $P_{\textrm{ult}}$ are measured in
$\textrm{GeV}/\textrm{fm}^{3}$ and $\textrm{MeV}/\textrm{fm}^{3}$,
respectively.

Based on Eq. (16) for $\varepsilon_{\textrm{ult}}$ and Eq. (17) for
$P_{\textrm{ult}}$, we find that the existence of a $2.08M_{\odot}$
($1.97M_{\odot}$) NS leads to $\varepsilon_{\textrm{c}}\lesssim
1.32\,\textrm{GeV}/\textrm{fm}^{3}$ ($1.47\,\textrm{GeV}/\textrm{fm}^{3}$) and
$P_{\textrm{c}}\lesssim 494\,\textrm{MeV}/\textrm{fm}^{3}$
($548\,\textrm{MeV}/\textrm{fm}^{3}$), respectively. It is necessary to point
out that the upper limit $P_{\textrm{c}}\lesssim
548\,\textrm{MeV}/\textrm{fm}^{3}$ is quite consistent with $P_{\max}\lesssim
570\,\textrm{MeV}/\textrm{fm}^{3}$ (shown as cyan diamond with error bars in
the inset of FIG. 9) obtained from Ref. [146] where the maximum masses are
assumed to be larger than $1.97M_{\odot}$ [185]. The general
$M_{\textrm{NS}}^{\max}$-dependence of $\varepsilon_{\textrm{ult}}$ and that
of $P_{\textrm{ult}}$ are shown by the magenta lines (together with the grey
band). Also shown are the predictions on $\varepsilon_{\textrm{ult}}$ and
$P_{\textrm{ult}}$ (black dashed lines with the lavender band) from Ref.
[131], i.e., $\varepsilon_{\textrm{ult}}\approx
7.62[M_{\odot}/M_{\textrm{NS}}^{\max}]^{2}\,\textrm{GeV}/\textrm{fm}^{3}$ and
$P_{\textrm{ult}}\approx
5.12[M_{\odot}/M_{\textrm{NS}}^{\max}]^{2}\,\textrm{GeV}/\textrm{fm}^{3}$
(cited as “LP, 2011”), respectively. Considering
$M_{\textrm{NS}}^{\max}/M_{\odot}\approx 2.08$, these relations lead to the
estimates $\varepsilon_{\textrm{ult}}\lesssim
1.76\,\textrm{GeV}/\textrm{fm}^{3}$ and $P_{\textrm{ult}}\lesssim
1.18\,\textrm{GeV}/\textrm{fm}^{3}$, which are about 33% and 139% larger than
those from Eq. (16) and Eq. (17), respectively. In addition, we find the
constraint $\varepsilon_{\textrm{c}}\lesssim
1.41\,\textrm{GeV}/\textrm{fm}^{3}$ (by using
$M_{\textrm{NS}}^{\max}/M_{\odot}\approx 2.01$ [186] in Eq. (16)) is close to
$\varepsilon_{\textrm{c}}\lesssim
1.18_{-0.17}^{+0.17}\,\textrm{GeV}/\textrm{fm}^{3}$ from recent ab-initio QCD
calculations [187], which also assumed the NS masses are greater than
$2.01M_{\odot}$ as shown by the blue hexagon. Similarly, we have the limit
$\varepsilon_{\textrm{c}}\lesssim 1.42\,\textrm{GeV}/\textrm{fm}^{3}$ (using
$M_{\textrm{NS}}^{\max}/M_{\odot}\approx 2$ in Eq. (16)), which is close to
the limit $1.46\,\textrm{GeV}/\textrm{fm}^{3}$ of Ref. [154] where the
algorithm rejects masses $\leq 2M_{\odot}$ (red circle). Furthermore, the
maximum central energy density $\lesssim
1.16_{-0.10}^{+0.11}\,\textrm{GeV}/\textrm{fm}^{3}$ from Ref. [188] is also
shown using the green diamond, in which the maximum mass
$M_{\textrm{NS}}^{\max}$ is greater than about 2.2$M_{\odot}$. Using our
estimate of Eq. (16), for such $M_{\textrm{NS}}^{\max}\approx 2.2M_{\odot}$
the central energy density is found to be about
$1.19\,\textrm{GeV}/\textrm{fm}^{3}$. We summarize in TAB. 1 these constraints
on $\varepsilon_{\textrm{ult}}$ and $P_{\textrm{ult}}$. Finally, by
considering the (black widow) PSR J0952-0607 [118] with its mass about
2.35$M_{\odot}$, the ultimate energy density $\varepsilon_{\textrm{ult}}$ and
pressure $P_{\textrm{ult}}$ are estimated to be about
$1.05\,\textrm{GeV}/\textrm{fm}^{3}$ and $392\,\textrm{MeV}/\textrm{fm}^{3}$
(lightblue bands), respectively.

##  VII Continuous Crossover of Sound Speed in Neutron Star Cores

Up to now, we have investigated a few consequences of the $s_{\textrm{c}}^{2}$
of Eq. (1) holding only at NS centers. In this section, we study the possible
occurrence of continuous crossover signaled by a smooth reduction of $s^{2}$
at high densities in NS cores (as sketched in panel (b) of FIG. 1), by
generalizing $s_{\textrm{c}}^{2}$ of Eq. (1) to finite distances away from NS
centers.

Figure 10: Sketch for two types of dependence of certain physical quantities
on energy density.

In this and next sections, we shall encounter two types of dependence of a
certain quantity $\mathcal{P}$ on the energy density, one is the dependence of
(central)
$\mathcal{P}_{\textrm{c}}\equiv\mathcal{P}(\varepsilon_{\textrm{c}})$ on
$\varepsilon_{\textrm{c}}$ (panel (a) of FIG. 10) or equivalently on
$\widehat{P}_{\textrm{c}}$, e.g., FIG. 11, FIG. 15, etc. The other is the
dependence of $\mathcal{P}$ on the reduced energy density
$\widehat{\varepsilon}=\varepsilon/\varepsilon_{\textrm{c}}$, e.g., Eq.(41)
and Eq. (VII.2), etc., see panel (b) of FIG. 10. The latter actually
encapsulates the radial-dependence of the quantity, i.e., a finite
$\widehat{r}$ corresponds to a finite $\widehat{\varepsilon}<1$. We may write
explicitly the subscripts to avoid potential confusions.

###  VII.1 1st-order in $\widehat{\varepsilon}-1$: Sign of
$s_{\textmd{c}}^{2}-s^{2}(\widehat{\varepsilon})$ for
$\widehat{\varepsilon}<1$

One can straightforwardly obtain the SSS at distance $\widehat{r}$ from the
center $\widehat{r}=0$ by considering the perturbative expressions for the
reduced pressure
$\widehat{P}\approx\widehat{P}_{\textrm{c}}+b_{2}\widehat{r}^{2}+b_{4}\widehat{r}^{4}+\cdots$
and the reduced energy density $\widehat{\varepsilon}\approx
1+a_{2}\widehat{r}^{2}+a_{4}\widehat{r}^{4}+\cdots$ [115],

$s^{2}(\widehat{r})=\frac{\mathrm{d}\widehat{P}}{\mathrm{d}\widehat{\varepsilon}}=\frac{\mathrm{d}\widehat{P}}{\mathrm{d}\widehat{r}}\frac{\mathrm{d}\widehat{r}}{\mathrm{d}\widehat{\varepsilon}}=\frac{b_{2}+2b_{4}\widehat{r}^{2}+\cdots}{a_{2}+2a_{4}\widehat{r}^{2}+\cdots}.$
(18)

Taking $\widehat{r}=0$ in Eq. (18) gives us $s_{\textrm{c}}^{2}=b_{2}/a_{2}$
and therefore $a_{2}=b_{2}/s_{\textrm{c}}^{2}<0$ since
$b_{2}=-6^{-1}(1+3\widehat{P}_{\textrm{c}}^{2}+4\widehat{P}_{\textrm{c}})<0$
[115], i.e., the energy density at the center is larger than its surroundings.
However, it does not necessarily mean that the center is stiffer than its
surroundings in terms of the SSS. We can expand Eq. (18) as an even power
series of $\widehat{r}$ (which becomes accurate as $\widehat{r}\to 0$),

$\displaystyle s^{2}(\widehat{r})\approx$ $\displaystyle
s_{\textrm{c}}^{2}\left[1+\frac{2}{b_{2}}\left(b_{4}-s_{\textrm{c}}^{2}a_{4}\right)\widehat{r}^{2}\right]+\mathcal{O}\left(\widehat{r}^{4}\right).$
(19)

The expression for $a_{4}$ could be obtained by using the expansion of
$\widehat{P}$ over $\widehat{\varepsilon}$. Specifically, by collecting the
terms in the expansion of
$\widehat{P}=\sum_{k=1}d_{k}\widehat{\varepsilon}^{k}$ in front of
$\widehat{r}^{4}$ and comparing it with $b_{4}\widehat{r}^{4}$ in the
expansion of $\widehat{P}$ over the coefficients $\\{b_{k}\\}$’s, one finds
that,

$a_{4}=\frac{1}{s_{\textrm{c}}^{2}}\left(b_{4}-a_{2}^{2}\sum_{k=1}^{K}\frac{k(k-1)}{2}d_{k}\right),$
(20)

where the following general sum rules are used,

$\displaystyle\widehat{P}_{\textrm{c}}=$
$\displaystyle\sum_{k=1}^{K}d_{k},~{}~{}s_{\textrm{c}}^{2}=\left.\frac{\mathrm{d}\widehat{P}}{\mathrm{d}\widehat{\varepsilon}}\right|_{\widehat{\varepsilon}=\widehat{\varepsilon}_{\textrm{c}}=1}=\sum_{k=1}^{K}kd_{k}.$
(21)

Here $K$ is the effective truncation order of the expansions. Putting the
expression of $a_{4}$ of Eq. (20) into Eq. (19) leads to the final expression
for $s^{2}$ as,

$\displaystyle s^{2}(\widehat{r})\approx
s_{\textrm{c}}^{2}+2a_{2}D\widehat{r}^{2},~{}~{}D=\sum_{k=1}^{K}\frac{k(k-1)}{2}d_{k}.$
(22)

Furthermore, since we have for the $\widehat{\varepsilon}$ away from the
center perturbatively as $\widehat{\varepsilon}(\widehat{r})\approx
1+a_{2}\widehat{r}^{2}$, the dependence of $s^{2}$ on $\widehat{r}^{2}$ could
be transformed into that on $\widehat{\varepsilon}$ as,

$s^{2}(\widehat{\varepsilon})\approx
s_{\textrm{c}}^{2}+2D\left(\widehat{\varepsilon}-1\right).$ (23)

The above approximation (23) is expected to be valid only near the NS centers
(nevertheless, the sign of $D$ is sufficient for our purpose to deduce a
relative relation between $s^{2}(\widehat{\varepsilon}$) and
$s_{\textrm{c}}^{2}$). According to Eq. (23), we then have

$\displaystyle D<0\;\leftrightarrow$
$\displaystyle\;s_{\textrm{c}}^{2}<s^{2}(\widehat{\varepsilon})$
$\displaystyle\leftrightarrow$ $\displaystyle\;\textrm{``reduction of
$s^{2}$\;toward\;NS\;centers}\mbox{''},$ (24)

since $\widehat{\varepsilon}-1<0$, see panel (b) of FIG. 1. In both Eq. (22)
and Eq. (23) the zeroth-order term is the SSS at center, what we investigate
is the possible crossover at zero temperature (cold NSs) and high densities
$\gtrsim 5\rho_{\textrm{sat}}$ [108, 109, 110, 111], with a peaked behavior in
$s^{2}$ (see next subsection).

Besides (21), the stability/causality condition for any
$\widehat{\varepsilon}$, i.e., $0\leq s^{2}\leq 1$, can be written as,

$0\leq
d_{1}+2d_{2}\widehat{\varepsilon}+3d_{3}\widehat{\varepsilon}^{2}+\cdots\leq
1,$ (25)

which guarantees the pressure $\widehat{P}$ never becomes negative. The
analyses given up to now are general and not limited to the maximum-mass
configuration $M_{\textrm{NS}}^{\max}$.

For the maximum-mass configuration $M_{\textrm{NS}}^{\max}$ on the NS M-R
curve, the inequality
$\mathrm{d}^{2}M_{\textrm{NS}}^{\max}/\mathrm{d}\varepsilon_{\textrm{c}}^{2}<0$
gives

$\displaystyle\left.\frac{\mathrm{d}s_{\textrm{c}}^{2}}{\mathrm{d}\widehat{P}_{\textrm{c}}}\right|_{M_{\textrm{NS}}^{\max}}<\sigma_{\textrm{c}}^{2}\equiv\frac{\mathrm{d}}{\mathrm{d}\widehat{P}_{\textrm{c}}}s_{\textrm{c}}^{2}=\frac{2}{3}\frac{9\widehat{P}_{\textrm{c}}^{4}-3\widehat{P}_{\textrm{c}}^{2}+4\widehat{P}_{\textrm{c}}+2}{(1-3\widehat{P}_{\textrm{c}}^{2})^{2}},$
(26)

here $s_{\textrm{c}}^{2}$ is given by Eq. (1). Inequality (26) implies, e.g.,
for PSR J0740+6620 that
$\mathrm{d}s_{\textrm{c}}^{2}/\mathrm{d}\widehat{P}_{\textrm{c}}\lesssim 2.74$
using $\widehat{P}_{\textrm{c}}\approx 0.24$ [115], which means if
$\widehat{P}_{\textrm{c}}$ increases by about 0.1, the increasing of
$s_{\textrm{c}}^{2}$ should be smaller than 0.274. Moreover, we obtain using
the relation
$\mathrm{d}s^{2}/\mathrm{d}\widehat{P}=\mathrm{d}s^{2}/\mathrm{d}\widehat{\varepsilon}\cdot\mathrm{d}\widehat{\varepsilon}/\mathrm{d}\widehat{P}=\mathrm{d}^{2}\widehat{P}/\mathrm{d}\widehat{\varepsilon}^{2}\cdot\mathrm{d}\widehat{\varepsilon}/\mathrm{d}\widehat{P}$
an equivalent form of (26),

$\left.\frac{\mathrm{d}^{2}\widehat{P}}{\mathrm{d}\widehat{\varepsilon}^{2}}\right|_{\widehat{\varepsilon}=\widehat{\varepsilon}_{\textrm{c}}=1}=\sum_{k=1}^{K}k(k-1)d_{k}=2D<\sigma_{\textrm{c}}^{2}s_{\textrm{c}}^{2}.$
(27)

The coefficients $d_{1}$ and $d_{2}$ (via conditions of (21)) and the
correction $D$ are given respectively by,

$\displaystyle d_{1}=$ $\displaystyle
2\widehat{P}_{\textrm{c}}-s_{\textrm{c}}^{2}+\sum_{k=3}^{K}(k-2)d_{k},$ (28)
$\displaystyle d_{2}=$
$\displaystyle-\widehat{P}_{\textrm{c}}+s_{\textrm{c}}^{2}-\sum_{k=3}^{K}(k-1)d_{k},$
(29) $\displaystyle D=$
$\displaystyle-3\widehat{P}_{\textrm{c}}+2s_{\textrm{c}}^{2}+\sum_{k=3}^{K}\frac{(k-2)(k-3)}{2}d_{k}.$
(30)

The first two terms in $D$ are deterministic while the last one has certain
randomness (reflecting the uncertainties of the dense matter EOS).

As an example, we use the above scheme to determine the correction coefficient
$D$ for $K=5$. The sum rules of (21) (or equivalently Eqs. (28) and (29)) can
be used to express the coefficients $d_{1}$ and $d_{2}$ as
$d_{1}=2\widehat{P}_{\textrm{c}}-s_{\textrm{c}}^{2}+d_{3}+2d_{4}+3d_{5}$ and
$d_{2}=-\widehat{P}_{\textrm{c}}+s_{\textrm{c}}^{2}-2d_{3}-3d_{4}-4d_{5}$,
respectively. Then, the vanishing of $s^{2}$ when $\widehat{\varepsilon}=0$ or
equivalently $d_{1}=0$ enables us to write
$d_{3}=-2\widehat{P}_{\textrm{c}}+s_{\textrm{c}}^{2}-2d_{4}-3d_{5}$.
Therefore, according to the basic definition of $D$ in Eq. (22), we obtain
$D=d_{2}+3d_{3}+6d_{4}+10d_{5}=d_{4}+3d_{5}+2s_{\textrm{c}}^{2}-3\widehat{P}_{\textrm{c}}$,
which is just the content of Eq. (30). See Appendix A for a proof of Eq. (30)
for general $K$.

Figure 11: The probability of $D<0$ for $K=3$ with $K$ the truncation order of
expansion
$\widehat{P}=\sum_{k=1}^{K}d_{k}\widehat{\varepsilon}^{k}=d_{1}\widehat{\varepsilon}+d_{2}\widehat{\varepsilon}^{2}+d_{3}\widehat{\varepsilon}^{3}$.

Figure 12: The same as FIG. 11 but for more cases of the truncation order $K$.
Dashed vertical blue line in each panel is for
$\widehat{P}_{\textrm{c}}\approx 0.24$ (PSR J0740+6620) extracted in Ref.
[115], and the black horizontal dashed line is plotted at
$\textrm{prob}(D<0)=50\%$.

The probability of $D<0$ could be estimated by combining the condition (25) of
$0\leq s^{2}\leq 1$, inequality (27) together with the general sum rules of
(21), i.e.,

$\displaystyle\textrm{prob}(D<0)\approx\frac{\\#[0\leq s^{2}\leq
1\;\textrm{and}\;2D<\sigma_{\textrm{c}}^{2}s_{\textrm{c}}^{2}\;\textrm{and}\;D<0]}{\\#[0\leq
s^{2}\leq 1\;\textrm{and}\;2D<\sigma_{\textrm{c}}^{2}s_{\textrm{c}}^{2}]},$
(31)

for the maximum-mass configuration $M_{\textrm{NS}}^{\max}$ of the NS M-R
curve, here $\\#$ counts for the number of samples of [$\dots$]. In the
following, we uniformly sample the coefficients $d_{k}$ for $k\geq 4$ within
certain ranges fulfilling the requirements $0\leq s^{2}\leq 1$ and
$2D<\sigma_{\textrm{c}}^{2}s_{\textrm{c}}^{2}$ and count the events of $D<0$
to estimate the probability of $D<0$. The value of $d_{2}$ could be obtained
correspondingly using Eq. (29) and $d_{3}$ via Eq. (28), i.e., the sum rules
of (21). For a certain truncation order $K$, different empirical ranges for
the parameters $d_{k}$ are adopted to boost the sampling efficiency.

The simplest situation is $K=3$, i.e., $\widehat{P}\approx
d_{1}\widehat{\varepsilon}+d_{2}\widehat{\varepsilon}^{2}+d_{3}\widehat{\varepsilon}^{3}\approx
d_{2}\widehat{\varepsilon}^{2}+d_{3}\widehat{\varepsilon}^{3}$.
Correspondingly, we shall obtain
$d_{2}=-s_{\textrm{c}}^{2}+3\widehat{P}_{\textrm{c}}$ and
$d_{3}=s_{\textrm{c}}^{2}-2\widehat{P}_{\textrm{c}}$, and the condition $D<0$
(see Eq. (30)) is now equivalently to (with the approximation holding for
small $\widehat{P}_{\textrm{c}}$ as we have $s_{\textrm{c}}^{2}\approx
4\widehat{P}_{\textrm{c}}/3$),

$u_{\textrm{c}}\equiv
2s_{\textrm{c}}^{2}-3\widehat{P}_{\textrm{c}}\approx-{\widehat{P}_{\textrm{c}}}/{3}<0.$
(32)

Considering Eq. (1), we then obtain $\widehat{P}_{\textrm{c}}\lesssim 0.105$
in order to make $D<0$ (see the inset of FIG. 11 where the diamond is at
$\widehat{P}_{\textrm{c}}\approx 0.105$). For $\widehat{P}_{\textrm{c}}\gtrsim
0.105$, there is no space for the occurrence of $D<0$, and the transition on
the probability of $D<0$ from $\widehat{P}_{\textrm{c}}\lesssim 0.105$ to
$\widehat{P}_{\textrm{c}}\gtrsim 0.105$ is sharp (FIG. 11). Equivalently, this
means $s_{\textrm{c}}^{2}$ is definitely larger than its surroundings for
$\widehat{P}_{\textrm{c}}\gtrsim 0.105$. Due to the importance of $D$ for our
analysis, we explain/illustrate in Appendix B why the $D$ tends to be negative
for small $\widehat{P}_{\textrm{c}}$ using $K=4$.

The above example with $K=3$ involving the cubic polynomial of $\widehat{P}$
over $\widehat{\varepsilon}$ is a very special case for demonstration. In
order to analyze the probability of $D<0$ more generally, we study more cases
of varying $K$ as shown in FIG. 12 (each panel simulates $10^{8}$ uniform
samples). It is seen that as the truncation order $K$ increases, the
probability of $D<0$ eventually stabilizes (see panels (a) to (h) of FIG. 12
where $K$ varies from 4 to 40). One finds from the curves that for
$\widehat{P}_{\textrm{c}}\lesssim 0.3$ the probability of $D<0$ is larger than
about 50%, indicating that the continuous crossover probably occurs near the
center. For PSR J0740+6620 with its mass $\approx 2.08\pm 0.07M_{\odot}$ [69,
70, 71, 72] which is close to the theoretical predicted maximum NS mass about
$2.01\mbox{$\sim$}2.16M_{\odot}$ [189] (Ref. [190] showed that it should be
smaller than about $2.17M_{\odot}$ while Ref. [191] predicted
$M_{\textrm{NS}}^{\max}\lesssim 2.16\mbox{$\sim$}2.28M_{\odot}$ via general
relativistic magnetohydrodynamics simulations) as well as
$\widehat{P}_{\textrm{c}}\approx 0.24$ and $s_{\textrm{c}}^{2}\approx 0.45$
[115], this probability is found to be larger than about 63%. As
$\widehat{P}_{\textrm{c}}$ increases even further, the probability of $D<0$
decreases. In the limiting case of $s_{\textrm{c}}^{2}\to 1$ (or equivalently
$\widehat{P}_{\textrm{c}}\to 0.374$), this probability is extremely small.

The physical reason for these behaviors could be traced back to the
(deterministic) term $2s_{\textrm{c}}^{2}-3\widehat{P}_{\textrm{c}}$ in $D$ of
Eq. (30) which is negative for small $\widehat{P}_{\textrm{c}}$ (see Eq.
(32)), though the summation term of Eq. (30) maybe either positive or negative
(EOS uncertainties). It indicates that $\widehat{P}_{\textrm{c}}$ is the
relevant quantity for the onset of the continuous crossover
($s^{2}(\widehat{\varepsilon})>s_{\textrm{c}}^{2}$) at zero temperature and
high densities (i.e., in the cores of cold NSs) [108, 109, 110, 111]: It is
likely to happen for large ($\widehat{P}_{\textrm{c}}$ is small and the core
is denser) and massive NSs (so they are near the maximum-mass configuration on
the M-R curve). At the Newtonian limit $s_{\textrm{c}}\approx
4\widehat{P}_{\textrm{c}}/3$, since
$u_{\textrm{c}}=-\widehat{P}_{\textrm{c}}/3$ the coefficient $D$ tends always
to be negative (e.g., see relevant discussion/illustration given in Appendix B
and FIG. b3). In this case, we expect that the probability of
$s^{2}(\widehat{\varepsilon})>s_{\textrm{c}}^{2}$ is large. Therefore the
general-relativistic effects reduce the probability for
$s^{2}(\widehat{\varepsilon})>s_{\textrm{c}}^{2}$.

Figure 13: Averages of $2a_{2}D/s_{\textrm{c}}^{2}$ and
$2a_{2}D_{-}/s_{\textrm{c}}^{2}$ as functions of $\widehat{P}_{\textrm{c}}$
(panel (a)) and coefficients $D$ and $J$ (panel (b)) defined in Eq. (22) and
Eq. (38), respectively, the vertical dashed line is at
$\widehat{P}_{\textrm{c}}=0.24$.

According to Eq. (22), the deviation of SSS from its central value is
proportional to $2a_{2}D\widehat{r}^{2}$. To evaluate the relative changing
rate of SSS from its central value, we show in panel (a) of FIG. 13 the
$\langle 2a_{2}D/s_{\textrm{c}}^{2}\rangle$ defined as

$\left\langle{2a_{2}D}/{s_{\textrm{c}}^{2}}\right\rangle=\left({2a_{2}}/{s_{\textrm{c}}^{2}}\right)\sum_{k=\pm}\textrm{prob}(D_{k})D_{k},$
(33)

here $D_{k}=D_{\pm}$ denotes the value of $D>0$ or $D<0$. For a comparison,
the average value of $2a_{2}D_{-}/s_{\textrm{c}}^{2}$ is also shown (blue
dashed line). Considering $\widehat{P}_{\textrm{c}}\approx 0.24$ [115] for PSR
J0740+6620 [69, 70, 71, 72], we obtain $\langle
2a_{2}D/s_{\textrm{c}}^{2}\rangle\approx 1.6$, i.e.,
$s^{2}(\widehat{r})/s_{\textrm{c}}^{2}\approx 1+1.6\widehat{r}^{2}$ for small
$\widehat{r}$. Transforming it back using Eq. (23) gives us
$s^{2}(\widehat{\varepsilon})/s_{\textrm{c}}^{2}\approx
1+2(1-\widehat{\varepsilon})\approx 3-2\widehat{\varepsilon}$, here
$\widehat{\varepsilon}$ being very close to 1 is assumed. Equivalently, we can
rewrite the radial-dependence of the SSS by using Eq. (6) approximately as (by
neglecting the intercept $0.64$),

${s^{2}}/{s_{\textrm{c}}^{2}}\approx
1+\frac{9.4\widehat{P}_{\textrm{c}}\langle
a_{2}D/s_{\textrm{c}}^{2}\rangle}{1+3\widehat{P}_{\textrm{c}}^{2}+4\widehat{P}_{\textrm{c}}}\left(\frac{r}{R_{\max}}\right)^{2},$
(34)

here $r=\widehat{r}Q$ is the physical radius and $R_{\max}$ is measured in km.
Taking $\widehat{P}_{\textrm{c}}\approx 0.24$, we obtain
$s^{2}/s_{\textrm{c}}^{2}\approx 1+0.85(r/R_{\max})^{2}$.

The consistency between Eq. (19) and Eq. (22) implies
$b_{4}-s_{\textrm{c}}^{2}a_{4}<0$ if $D<0$ due to the negativeness of $a_{2}$
and $b_{2}$. Consequently, we obtain $a_{4}>{b_{4}}/{s_{\textrm{c}}^{2}}>0$
since $b_{4}>0$ [115]. As a rough (order of magnitude) estimate for $a_{4}$,
e.g., we take $\widehat{P}_{\textrm{c}}\gtrsim 0.16$ (for massive NSs),
therefore $a_{4}>b_{4}/s_{\textrm{c}}^{2}\gtrsim 0.86$, i.e.,
$a_{4}\sim\mathcal{O}(1)$ (see Appendix D for another order of magnitude
estimate for $a_{4}$). Since all of our perturbative expansions including
$\widehat{\varepsilon}$ and $\widehat{\rho}=\rho/\rho_{\textrm{c}}$ (see Eq.
(53)) over $\widehat{r}$ are dimensionless, it is reasonable to expect the
magnitude of the expansion coefficients is natural, i.e., it is
$\sim\mathcal{O}(1)$. On the other hand, if $D>0$, then the consistency
between Eq. (19) and Eq. (22) implies $a_{4}<b_{4}/s_{\textrm{c}}^{2}$.
Combining these two cases leads to

$D\left({b_{4}}/{s_{\textrm{c}}^{2}}-a_{4}\right)>0,$ (35)

i.e., $D$ and $b_{4}/s_{\textrm{c}}^{2}-a_{4}$ have the same sign. The
inequality (35) is a general constraint on the expansion coefficients.

###  VII.2 2nd-order in $\widehat{\varepsilon}-1$: Estimates on the Peak of
$s^{2}$

The $s^{2}$ approaches zero at the core-crust transition density below which
the dynamical fluctuations in uniform matter will be unstable indicating the
onset of forming clusters in the crust. It implies that there should exist a
peak in $s^{2}$ somewhere along the path away from the center if $D<0$. In
fact, Refs. [151, 152] predicted that the peak of $s^{2}$ for NSs with masses
$M_{\textrm{NS}}^{\max}\gtrsim 2M_{\odot}$ is located at about 2/3 of its
radius. As $\varepsilon$ is a decreasing function of $\widehat{r}$, the peak
of $s^{2}$ is located at a energy density
$\varepsilon_{\textrm{pk}}<\varepsilon_{\textrm{c}}$, e.g., Ref. [192]
predicted that $\varepsilon_{\textrm{pk}}\approx
565\,\textrm{MeV}/\textrm{fm}^{3}$ for massive NSs. Take our over-simplified
EOS with truncation order $K=3$ (FIG. 11) as an example, the peak
$\widehat{\varepsilon}_{\textrm{pk}}$ of $s^{2}(\widehat{\varepsilon})$ can be
obtained as $\widehat{\varepsilon}_{\textrm{pk}}\approx
5/6+\widehat{P}_{\textrm{c}}+7\widehat{P}_{\textrm{c}}^{2}/2+13\widehat{P}_{\textrm{c}}^{3}\gtrsim
5/6$ and $s^{2}(\widehat{\varepsilon}_{\textrm{pk}})\approx
25\widehat{P}_{\textrm{c}}/18+5\widehat{P}_{\textrm{c}}^{2}/9+17\widehat{P}_{\textrm{c}}^{3}/6$
(with $\widehat{P}_{\textrm{c}}$ being small $\lesssim 0.105$) using
$s^{2}\approx
2(3\widehat{P}_{\textrm{c}}-s_{\textrm{c}}^{2})\widehat{\varepsilon}+3(s_{\textrm{c}}^{2}-2\widehat{P}_{\textrm{c}})\widehat{\varepsilon}^{3}$.
In addition, the peak of $s^{2}$ could not be estimated through Eq. (23) since
it is linear in $\widehat{\varepsilon}$ and the higher-order terms, e.g.,
$a_{4}$, $b_{4}$ and $a_{6}$, etc., are necessary for such estimate.

We now use our perturbative expansions as a tool and include these higher-
order contributions to explore whether there would exist a plateau for
$\widehat{P}$ for a finite range of $\widehat{\varepsilon}$ (relevant for
sharp PTs) and estimate the peak of $s^{2}$ (near NS centers) if it exists.
When considering terms including order $b_{6}$-coefficient in the expansion of
$\widehat{P}$ and the $a_{6}$-term in the expansion of
$\widehat{\varepsilon}$, namely
$\widehat{P}\approx\widehat{P}_{\textrm{c}}+b_{2}\widehat{r}^{2}+b_{4}\widehat{r}^{4}+b_{6}\widehat{r}^{6}$
and $\widehat{\varepsilon}\approx
1+a_{2}\widehat{r}^{2}+a_{4}\widehat{r}^{4}+a_{6}\widehat{r}^{6}$, we derive
the SSS like Eq. (19),

$\displaystyle s^{2}/s_{\textrm{c}}^{2}\approx$ $\displaystyle
1+\frac{2}{b_{2}}\left(b_{4}-s_{\textrm{c}}^{2}a_{4}\right)\widehat{r}^{2}$
$\displaystyle+\frac{3}{b_{2}}\left[\left(b_{6}-s_{\textrm{c}}^{2}a_{6}\right)-\frac{4}{3}\frac{a_{4}}{a_{2}}\left(b_{4}-s_{\textrm{c}}^{2}a_{4}\right)\right]\widehat{r}^{4}.$
(36)

The coefficient $a_{6}$ could be eliminated as $a_{4}$ in Eq. (19): Expanding
the reduced pressure
$\widehat{P}=\sum_{k=1}^{K}d_{k}\widehat{\varepsilon}^{k}$ with
$\widehat{\varepsilon}=1+\sum_{k=1}^{K}a_{k}\widehat{r}^{k}$ and comparing the
coefficient of $\widehat{r}^{6}$ in the expansion
$\widehat{P}=\widehat{P}_{\textrm{c}}+\sum_{k=1}^{K}b_{k}\widehat{r}^{k}$. The
result is,

$b_{6}-s_{\textrm{c}}^{2}a_{6}=a_{2}^{3}J+2a_{2}a_{4}D,$ (37)

where $J$ is a quantity constructed from $d_{k}$’s defined in the expansion
$\widehat{P}=\sum_{k=1}^{K}d_{k}\widehat{\varepsilon}^{k}$,

$J=\sum_{k=1}^{K}\frac{k(k-1)(k-2)}{6}d_{k}=d_{3}+4d_{4}+10d_{5}+\cdots.$ (38)

The coefficient $J$ has certain randomness like the coefficient $D$
(characterizing the uncertainties of the dense matter EOS). Using this $J$, we
can rewrite the radial variation of $s^{2}$ as,

$s^{2}(\widehat{r})\approx
s_{\textrm{c}}^{2}+2a_{2}D\widehat{r}^{2}+\left(3a_{2}^{2}J+2a_{4}D\right)\widehat{r}^{4}.$
(39)

As discussed above, the peak of $s^{2}$ exists in the situation of $D<0$ and
where $a_{4}>b_{4}/s_{\textrm{c}}^{2}>0$ (see inequality (35)). Taking the
derivative $\mathrm{d}s^{2}/\mathrm{d}\widehat{r}$ and setting it to zero
locates the peak of $s^{2}$ to

$\widehat{r}_{\textrm{pk}}=\left(-\frac{a_{2}D}{3a_{2}^{2}J+2a_{4}D}\right)^{1/2}.$
(40)

Here $a_{2}D>0$, so the expression under square root is positive if
$J<-2a_{4}D/3a_{2}^{2}$ (with the latter being positive since $a_{4}D<0$). A
further quantitative estimate for $\widehat{r}_{\textrm{pk}}$ depends strongly
on the coefficients $a_{4}$, $D$ and $J$. Inverting
$\widehat{\varepsilon}\approx
1+a_{2}\widehat{r}^{2}+a_{4}\widehat{r}^{4}+a_{6}\widehat{r}^{6}$ gives us
$\widehat{r}^{2}\approx(\mu/a_{2})[1-(a_{4}/a_{2}^{2})\mu+(2a_{4}^{2}/a_{2}^{4}-a_{6}/a_{2}^{3})\mu^{2}]$,
where $\mu\equiv\widehat{\varepsilon}-1<0$. Putting this $\widehat{r}^{2}$
into the full expression for $s^{2}(\widehat{r})$ to order $\widehat{r}^{4}$
of Eq. (39) enables us to write the $s^{2}(\widehat{\varepsilon})=s^{2}(\mu)$
as,

$\displaystyle s^{2}(\mu)\approx$ $\displaystyle
s_{\textrm{c}}^{2}+2D\mu+3J\mu^{2}-\frac{2}{a_{2}^{2}}\left({3a_{4}}J+\frac{a_{6}}{a_{2}}D\right)\mu^{3},$
(41)

which generalizes (23). Thus we have to order $\mathcal{O}(\mu^{3})$,

$\displaystyle\widehat{P}(\mu)\approx$
$\displaystyle\widehat{P}_{\textrm{c}}+s_{\textrm{c}}^{2}\mu+\left.\frac{1}{2}\frac{\mathrm{d}s^{2}}{\mathrm{d}\mu}\right|_{\mu=0}\mu^{2}+\left.\frac{1}{6}\frac{\mathrm{d}^{2}s^{2}}{\mathrm{d}\mu^{2}}\right|_{\mu=0}\mu^{3}$
$\displaystyle\approx$
$\displaystyle\widehat{P}_{\textrm{c}}+s_{\textrm{c}}^{2}\mu+D\mu^{2}+J\mu^{3}.$
(42)

In fact, Eq. (VII.2) could be obtained straightforwardly via the expansion
$\widehat{P}=\sum_{k=1}^{K}d_{k}\widehat{\varepsilon}^{k}$. However, the
coefficient $J$ (Eq. (38)) now like the coefficient $D$ (Eq. (30)) is not
totally random, but is constrained through the general requirements $0\leq
s^{2}\leq 1$ and $\sigma_{\textrm{c}}^{2}s_{\textrm{c}}^{2}>2D$ used in the
algorithm of Eq. (31), i.e., the randomness of $J$ and $D$ is effectively
reduced by the physical constraints. Therefore, we can estimate $J$ in a
similar manner as for $D$. Taking $\widehat{P}_{\textrm{c}}\approx 0.24$ (for
PSR J0740+6620) we obtain an estimate $J\approx-2.7$ (panel (b) of FIG. 13).
Combining $s_{\textrm{c}}^{2}\approx 0.45$ and $D\approx-0.45$ (panel (b) of
FIG. 13) for $\widehat{P}_{\textrm{c}}\approx 0.24$, the
$\widehat{P}(\widehat{\varepsilon})$ is found to be an increasing function of
$\widehat{\varepsilon}$ for $\widehat{\varepsilon}\gtrsim 0.7$. It basically
excludes the possibility of having a plateau in $\widehat{P}$ at high
densities near the center. Including the last term in Eq. (41) may slightly
shift $\widehat{\varepsilon}\gtrsim 0.7$ to a higher value, here $a_{6}$ could
be estimated via Eq. (37) with the expression for coefficient $b_{6}$ given in
Ref. [115]. Moreover, the $D$-term and $J$-term contribute to $s^{2}(\mu)$
(and also $\widehat{P}(\mu)$) with different signs (since $\mu<0$) for
$\widehat{P}_{\textrm{c}}\lesssim 0.3$ (panel (b) of FIG. 13), explaining the
existence of a peak in $s^{2}(\mu)$. One finds the peak of $s^{2}(\mu)$ from
Eq. (41) (neglecting the last term) as,

$\mu_{\textrm{pk}}\approx-D/3J.$ (43)

We can see that the possible peak in the derivative part of $s^{2}$ is shifted
to a slightly lower value (see Eq. (50)).

As an illustration, we take $\varepsilon_{\textrm{c}}\approx
901\,\textrm{MeV}/\textrm{fm}^{3}$ (which is the central energy density for
PSR J0740+6620) [115] and $R\approx 12.39\,\textrm{km}$ [70] while let
$\widehat{P}_{\textrm{c}}$ be a free parameter to demonstrate the numerical
results. The radius length is obtained as $Q=(4\pi
G\varepsilon_{\textrm{c}})^{-1/2}\approx 8.7\,\textrm{km}$. Shown in FIG. 14
are our estimates for the peak in $s^{2}$ for either
$s^{2}(\widehat{\varepsilon})$ or $s^{2}(\widehat{r})$. Specifically,
$\widehat{\varepsilon}_{\textrm{pk}}=1+\mu_{\textrm{pk}}\approx 1-D/3J$ and
$\widehat{r}_{\textrm{pk}}$ is given by Eq. (40). The enhancement on $s^{2}$,
i.e., $\Delta
s^{2}=s^{2}(\widehat{r}_{\textrm{pk}})-s_{\textrm{c}}^{2}=-D^{2}/(3J+2a_{4}D/a_{2}^{2})$
from Eq. (39) is also shown, where $a_{4}\approx 1$ is adopted. Considering
that $|2a_{4}D/a_{2}^{2}|\ll|3J|$, we have $\Delta s^{2}\approx-D^{2}/3J$
using Eq. (43) (see the cyan dash-dotted line in FIG. 14 for
$15\cdot(-D^{2}/3J)=-5D^{2}/J$). The difference between $\Delta s^{2}$ using
Eq. (39) and Eq. (41) is that Eq. (39) is full to order $\widehat{r}^{4}$
while Eq. (41) is an approximation (even the $\mu^{3}$-term is included). It
is seen that as the $\widehat{P}_{\textrm{c}}$ increases, the
$\widehat{\varepsilon}_{\textrm{pk}}$ eventually approaches 1, implying that
the peak $\widehat{\varepsilon}_{\textrm{pk}}$ eventually moves to the center,
and it may even disappear if the $\widehat{P}_{\textrm{c}}$ increases further.
Actually for $\widehat{P}_{\textrm{c}}\gtrsim 0.3$ (see panel (b) of FIG. 13),
the $D$ changes from being negative to positive and therefore both terms
$2D\mu$ and $3J\mu^{2}$ in $s^{2}(\mu)$ are negative, i.e.,
$s_{\textrm{c}}^{2}$ is larger than its surroundings. Similar phenomena occur
for the $\Delta s^{2}$ and $\widehat{r}_{\textrm{pk}}Q/R$. For example, we
have for $\widehat{P}_{\textrm{c}}\to 0.374$ that $s_{\textrm{c}}^{2}\to 1$
(Eq. (1)) and therefore $\Delta s^{2}\to 0$ (no space under this limit for
$s^{2}$ to be enhanced when going outward from the center). In our
illustration, we find $(\Delta s^{2})_{\max}\approx 0.038$ occurring for
$\widehat{P}_{\textrm{c}}\approx 0.16$, therefore $(\Delta
s^{2})_{\max}/s_{\textrm{c}}^{2}\approx 14\%$ using $s_{\textrm{c}}^{2}\approx
0.26$. Similarly, we have $\Delta s^{2}/s_{\textrm{c}}^{2}\approx 5\%$ for
$\widehat{P}_{\textrm{c}}\approx 0.24$.

Figure 14: Estimate for the peak in $s^{2}$ either for
$s^{2}(\widehat{\varepsilon})$ or $s^{2}(\widehat{r})$, the enhancement on
$s^{2}$, i.e., $\Delta
s^{2}=s^{2}(\widehat{r}_{\textrm{pk}})-s_{\textrm{c}}^{2}=-D^{2}/(3J+2a_{4}D/a_{2}^{2})$
and its approximation $15\cdot(-D^{2}/3J)=-5D^{2}/J$ are shown. The
coefficient $D$ becomes averagely negative when
$\widehat{P}_{\textrm{c}}\gtrsim 0.3$ (see FIG. 13). In the simulation, the
effective truncation order of the expansions is set at $K=8$.

For situations with positive $D$ (thus $a_{2}D<0$), the SSS may firstly
decrease as finite $\widehat{r}$ develops and then increase further if
$3a_{2}^{2}J+2a_{4}D>0$. There exists then a valley
$\widehat{r}_{\textrm{va}}$ in $s^{2}$, which may still be estimated by using
Eq. (40). The self-consistent requirement between expansions (VII.2) and (39)
considering the coefficients of $\widehat{r}^{4}$ gives us (see the similar
inequality (35) from the preceding order),

$\left[\frac{4a_{4}}{3a_{2}}\left(\frac{b_{4}}{s_{\textrm{c}}^{2}}-a_{4}\right)-\left(\frac{b_{6}}{s_{\textrm{c}}^{2}}-a_{6}\right)\right]\left(J+\frac{2a_{4}}{3a_{2}^{2}}D\right)>0,$
(44)

where $b_{4}$ and $b_{6}$ were given in the Appendix of Ref. [115] and $a_{4}$
as well as $a_{6}$ are expected to be $\mathcal{O}(1)$.

The above results/estimates on the peak of $s^{2}$ as well as the qualitative
predictions for the $\widehat{P}(\widehat{\varepsilon})=\widehat{P}(\mu)$ are
expected to be (nearly) model independent, because we rely only on (without
using extra assumptions) the general conditions/requirements (including the
inequality $0\leq s^{2}\leq 1$ of (25), inequality (27) and the sum rules of
(21)) together with formula (1) and also bound the high-order coefficients
such as $a_{4}$ (Appendix D) and $a_{6}$, etc.

It is necessary to emphasize that due to its perturbative nature, however, the
above analysis on the $\widehat{\varepsilon}$-dependence of $s^{2}$ (or the
signature of the continuous crossover) is only effective and qualitatively
reasonable at finite distances $\widehat{r}$ being very close to NS centers
($\widehat{r}=0$). Whether there exist sharp PTs at some large distances away
from the centers (or at certain intermediate energy density $\varepsilon$)
could not be excluded immediately and should be in principle analyzed by
including more higher-order terms in the expansions, and the results may have
certain model dependence.

Furthermore, $\widehat{P}(\mu)<\widehat{P}_{\textrm{c}}$ is definite for
$\mu\approx 0^{-}$ considering Eq. (VII.2) which is fundamentally different
from $s^{2}(\mu)$ of Eq. (41), since $D$ still has sizable probabilities to be
positive even for small $\widehat{P}_{\textrm{c}}$ and thus
$s^{2}(\mu)<s_{\textrm{c}}^{2}$ (panel (c) of FIG. 1). We discuss in further
details in the next section on the possible bound on $P/\varepsilon$ in NS
cores.

##  VIII Upper Bound for $P/\varepsilon$ and Conformal Anomaly near Neutron
Star Centers

In this section, we discuss more on $P/\varepsilon$ in NS cores. Considering
the trace anomaly (or the conformality measure) $\Delta=1/3-P/\varepsilon$ and
compared with the analysis given in Ref. [119], our formula (1) predicts a
slightly different bound for $\Delta_{\textrm{c}}$ as
$\Delta_{\textrm{c}}\gtrsim-0.04$, and for $\widehat{P}_{\textrm{c}}\approx
0.24$ (PSR J0740+6620) we found $\Delta_{\textrm{c}}\approx 0.09$. Although
our $\Delta_{\textrm{c}}$ and that from TA analysis are close to each other,
the physical origins seem different. While the TA is based on properties of
pQCD, ours could be traced back to the structures encapsulated in the
(general-relativistic) TOV equations. Thus, whether there exists deep
connection between them needs further investigations.

Firstly, we mention a few similar evaluations of $\Delta$ in the recent
literature. An analysis in Ref. [152] using an agnostic EOS showed that
$\Delta$ is very close to zero for $M_{\textrm{NS}}^{\max}\gtrsim
2.18\mbox{$\sim$}2.35M_{\odot}$ and may be slightly negative for even more
massive NSs (e.g., $\Delta\gtrsim-0.02$ for $M_{\textrm{NS}}^{\max}\gtrsim
2.52M_{\odot}$). Moreover, in Ref. [192], the central minimum value of
$\Delta$ is found to be about $\Delta_{\min}\approx 0.04$ using the NICER data
together with the tidal deformability from GW170817. A value of
$\Delta_{\min}\approx-0.05$ was inferred considering additionally the second
component of GW190814 as a neutron star with mass about $2.59M_{\odot}$ [163]
using two hadronic EOS models [192]. See also the recent works on the similar
issue in Refs. [193, 194], e.g., Ref. [193] predicted that
$\Delta_{\textrm{c}}\gtrsim-0.046$ if $M_{\textrm{NS}}^{\max}\geq
2.2M_{\odot}$. Another analysis within the Bayesian framework considering the
state-of-the-art theoretical calculations showed that $\Delta\gtrsim-0.01$
(where $M_{\textrm{NS}}^{\max}\approx 2.27_{-0.11}^{+0.11}M_{\odot}$) and for
a $2M_{\odot}$ NS the polytropic index $\gamma\equiv\mathrm{d}\ln
P/\mathrm{d}\ln\varepsilon=s^{2}/(P/\varepsilon)$ is found to be about 2
[195]. The $\gamma$ parameter in NS centers, actually, can not be unity [119,
195] since $\Delta_{\textrm{c}}$ ($\widehat{P}_{\textrm{c}}$) and
$s_{\textrm{c}}^{2}$ could not approach zero (1/3) and 1/3 simultaneously, as
shown by Eq. (1). Physically, due to the nonlinear nature of the central EOS
in NSs, it is likely that $\gamma\neq 1$ (actually $\gamma=1$ if and only if
$P\propto\varepsilon$), see the relevant discussions in detail in Section IV
and FIG. 7. Using our formula (1), it is straightforward to find that
$4/3\leq\gamma_{\textrm{c}}\lesssim 2.67$ and considering
$\widehat{P}_{\textrm{c}}\approx 0.24_{-0.07}^{+0.05}$ for PSR J0740+6620
[115], we obtain $\gamma_{\textrm{c}}\approx 1.86_{-0.21}^{+0.22}$. Moreover,
for $\Delta_{\textrm{c}}=0$ the $\gamma_{\textrm{c}}$ index is $7/3\approx
2.33$. Furthermore, Ref. [195] also defined a quantity
$\Theta\equiv[\Delta^{2}+(P/\varepsilon-s^{2})^{2}]^{1/2}$ to measure the
conformality and found it should be $\lesssim 0.65$ for all NSs. The criterion
$\Theta\lesssim 0.2$ (besides $\gamma\lesssim 1.75$) was adopted [195] to
identify the nearly-conformality at a given density. Based on our formula (1),
we can find $0.19\lesssim\Theta_{\textrm{c}}\lesssim 0.63$ with the
minimum/maximum value obtained at $\widehat{P}_{\textrm{c}}\approx
0.18$/$\widehat{P}_{\textrm{c}}\approx 0.374$ and $\Theta_{\textrm{c}}\approx
0.22_{-0.03}^{+0.09}$ for $\widehat{P}_{\textrm{c}}\approx
0.24_{-0.07}^{+0.05}$ (PSR J0740+6620).

Quantity | Range | $\widehat{P}_{\textrm{c}}\approx 0.18$ | $0.24$ | $0.374$
---|---|---|---|---
$\Delta_{\textrm{c}}=1/3-\widehat{P}_{\textrm{c}}$ | $-0.041\lesssim\Delta_{\textrm{c}}\leq 1/3$ | $0.15$ | $0.09$ | $-0.041$
$s_{\textrm{c}}^{2}=\mathrm{d}{P}_{\textrm{c}}/\mathrm{d}\varepsilon_{\textrm{c}}$ | $0\leq s_{\textrm{c}}^{2}\leq 1$ | $0.31$ | $0.45$ | $1$
$\gamma_{\textrm{c}}=s_{\textrm{c}}^{2}/\widehat{P}_{\textrm{c}}$ | $4/3\leq\gamma_{\textrm{c}}\lesssim 2.67$ | $1.68$ | 1.86 | $2.67$
$t_{\textrm{c}}=\widehat{P}_{\textrm{c}}-s_{\textrm{c}}^{2}$ | $-0.63\lesssim t_{\textrm{c}}\leq 0$ | $-0.13$ | $-0.21$ | $-0.63$
$\Theta_{\textrm{c}}=(\Delta_{\textrm{c}}^{2}+t_{\textrm{c}}^{2})^{1/2}$ | $0.19\lesssim\Theta_{\textrm{c}}\lesssim 0.63$ | $0.19$ | 0.22 | $0.63$
$u_{\textrm{c}}=2s_{\textrm{c}}^{2}-3\widehat{P}_{\textrm{c}}$ | $-0.01\lesssim u_{\textrm{c}}\lesssim 0.88$ | $0.07$ | $0.17$ | $0.88$
TAB. 2: Ranges of several quantities relevant for the trace anomaly as defined
in the text, the values for them at three reference $\widehat{P}_{\textrm{c}}$
are also shown (last three columns),
$t_{\textrm{c}}=\mathrm{d}\Delta_{\textrm{c}}/\mathrm{d}\ln\varepsilon_{\textrm{c}}=\widehat{P}_{\textrm{c}}-s_{\textrm{c}}^{2}$
is the logarithmic derivative of $\Delta_{\textrm{c}}$, and
$u_{\textrm{c}}=2s_{\textrm{c}}^{2}-3\widehat{P}_{\textrm{c}}$ is the quantity
closely related to the possible crossover occurring near the NS centers (see
Eq. (30)).  Figure 15: Decomposition of $s_{\textrm{c}}^{2}$ into its non-
derivative part ($\widehat{P}_{\textrm{c}}$) and derivative part
($s_{\textrm{c}}^{2}-\widehat{P}_{\textrm{c}}$) for
$M_{\textrm{NS}}^{\max}/M_{\odot}=2$. Solid diamonds characterize the boundary
of $s_{\textrm{c}}^{2}\leq 1$ (left of the lavender band).

We summarize in TAB. 2 these values at three reference
$\widehat{P}_{\textrm{c}}$ (0.18, 0.24 and 0.374), where
$t_{\textrm{c}}\equiv\mathrm{d}\Delta_{\textrm{c}}/\mathrm{d}\ln\varepsilon_{\textrm{c}}=\widehat{P}_{\textrm{c}}-s_{\textrm{c}}^{2}$
as the logarithmic derivative of $\Delta_{\textrm{c}}$ [195] (with respect to
$\varepsilon_{\textrm{c}}$) is also given. We have $t_{\textrm{c}}\to 0$ if
$\gamma_{\textrm{c}}\to 1$ (the so-called conformal limit of matter),
therefore $-t_{\textrm{c}}$ characterizes the deviation from the conformal
limit. In Ref. [119], $-t=s^{2}-P/\varepsilon$ and $1/3-\Delta=P/\varepsilon$
are decomposed as the derivative and non-derivative parts of $s^{2}$,
respectively. Equivalently, one has $\Theta=(\Delta^{2}+t^{2})^{1/2}$. See
FIG. 15 for an example where $M_{\textrm{NS}}^{\max}/M_{\odot}=2$ is adopted
for the illustration (using Eq. (7) to solve for $\widehat{P}_{\textrm{c}}$
and $s_{\textrm{c}}^{2}$). According to Eq. (1), we find both the non-
derivative and derivative parts of $s_{\textrm{c}}^{2}$ are increasing
functions of $\varepsilon_{\textrm{c}}$, and take their maximum values at
$\widehat{P}_{\textrm{c}}\approx 0.374$. The solid diamonds in FIG. 15
characterize the boundary $\widehat{P}_{\textrm{c}}\lesssim 0.374$, i.e., the
curves above $\varepsilon_{\textrm{c}}\gtrsim
1.42\,\textrm{GeV}/\textrm{fm}^{3}$ (for $M_{\textrm{NS}}^{\max}/M_{\odot}=2$,
see TAB. 1) violate the causality condition $s_{\textrm{c}}^{2}\leq 1$. For
other values of $M_{\textrm{NS}}^{\max}$ the shapes of the curves in FIG. 15
are similar. Shortly after, we may show that there exists a peak in the
derivative part of $s^{2}$ when considering positions having finite distance
from the center (see Eq. (50)).

Figure 16: Dependence of $\Theta_{\textrm{c}}$ on the polytropic index
$\gamma_{\textrm{c}}$ for the maximum-mass configuration
$M_{\textrm{NS}}^{\max}$, where in the lightblue region two criteria
$\Theta_{\textrm{c}}\lesssim 0.2$ and $\gamma_{\textrm{c}}\lesssim 1.75$ hold
simultaneously. The constraints from Ref. [195] on $\gamma$ and $\Theta$ are
also shown (green solid lines associated with arrows) for comparison. The
inset amplifies the lightblue region.

In TAB. 2, the column with $\widehat{P}_{\textrm{c}}\approx 0.18$ is also
shown as the measure $\Theta_{\textrm{c}}$ takes its minimum at this
$\widehat{P}_{\textrm{c}}$. Interestingly, if one requires
$\Theta_{\textrm{c}}\approx 0.2$, then two solutions
$\widehat{P}_{\textrm{c}}\approx 0.21$ and $\widehat{P}_{\textrm{c}}\approx
0.16$ should be obtained, the corresponding $\gamma_{\textrm{c}}$ is found to
be about 1.75 and 1.62, respectively. We have approximately
$\Theta_{\textrm{c}}\approx
1/3-\eta+2\eta^{2}/3+13\eta^{3}/6+41\eta^{4}/8+\cdots$, where
$0.5\gtrsim\eta\equiv 1-4/3\gamma_{\textrm{c}}\geq 0$ acts as a small-
expansion quantity (since $2.67\gtrsim\gamma_{\textrm{c}}\geq 4/3$). These
numerical values indicate that empirically the near-conformality may be
possible in the core of very massive NSs (with masses $\gtrsim 2M_{\odot}$).
We notice that $\widehat{P}_{\textrm{c}}\approx 0.189$ [115] and
correspondingly $\Theta_{\textrm{c}}\approx 0.195$ and
$\gamma_{\textrm{c}}\approx 1.70$ are obtained if the radius about
$13.7\,\textrm{km}$ was adopted for PSR J0740+6620 [71]. However, the
likelihood for realizing conformality near the centers is small if the
criteria $\Theta_{\textrm{c}}\lesssim 0.2$ and $\gamma_{\textrm{c}}\lesssim
1.75$ [195] were adopted simultaneously, as shown in FIG. 16, where the
constraints $\gamma\lesssim 1.79$ and $\Theta\lesssim 0.24$ from Ref. [195]
are shown (green solid lines associated with arrows) for comparison. Moreover,
Eq. (55)$\sim$Eq. (59) given in the following together may extend this
conclusion to situations where $\widehat{\varepsilon}\lesssim 1$ (or
equivalently at some finite distances away from the center). Furthermore, FIG.
16 also supports the conclusion that the criterion
$\Theta_{\textrm{c}}\lesssim 0.2$ is more restrictive than
$\gamma_{\textrm{c}}\lesssim 1.75$ [195].

In FIG. 17, we show the $\Delta_{\textrm{c}}=1/3-\widehat{P}_{\textrm{c}}$ as
a function of $\varepsilon_{\textrm{c}}$ for different values of
$M_{\textrm{NS}}^{\max}$ (varying from 1.7$M_{\odot}$ to $2.4M_{\odot}$) using
the mass correlation (6). The parametrization for $\Delta$ suggested by Ref.
[119] is also shown (green dotted line), which is positive-definite by
construction. In addition, two model predictions for $\Delta$ are represented
by lavender [119] and cyan dashed [187] bands using machine-learning
algorithms and ab-initio QCD calculations, respectively. For relatively light
NSs, the condition $\Delta_{\textrm{c}}\geq 0$ generally holds since the
central energy density $\varepsilon_{\textrm{c}}$ is relatively low, e.g., for
$M_{\textrm{NS}}^{\max}/M_{\odot}=1.7$ the condition
$\varepsilon_{\textrm{c}}\lesssim 1.8\,\textrm{GeV}/\textrm{fm}^{3}$ is safely
satisfied. On the other hand, the $\varepsilon_{\textrm{c}}$ may exceed the
value set by $\Delta_{\textrm{c}}\approx 0$ as $M_{\textrm{NS}}^{\max}$
increases. Our constraint on the lower bound of $\Delta_{\textrm{c}}$ is quite
consistent with those from Refs. [119, 187]. The inset of FIG. 17 shows the
dependence of $\Delta_{\textrm{c}}$ on $s_{\textrm{c}}^{2}$. It is seen that
for $s_{\textrm{c}}^{2}\geq 7/9\approx 0.78$ (see Eq. (1)) the
$\Delta_{\textrm{c}}$ becomes negative. In short, FIG. 17 clearly explains why
the criterion $\Delta\geq 0$ tends to break down for massive NSs with
increasing $M_{\textrm{NS}}^{\max}$ [152, 192, 195, 196], instead of for light
NSs.

Figure 17: Conformality measure
$\Delta_{\textrm{c}}=1/3-\widehat{P}_{\textrm{c}}$ as a function of
$\varepsilon_{\textrm{c}}$ for different $M_{\textrm{NS}}^{\max}/M_{\odot}$
(from 1.7 to 2.4), where lightblue solid diamonds mark the points set by
$\widehat{P}_{\textrm{c}}\lesssim 0.374$ while the vertical black dashed lines
are for the vanishing points of $\Delta_{\textrm{c}}$ for each
$M_{\textrm{NS}}^{\max}$. Two model calculations (lavender/cyan) [119, 187]
and a parametrization for $\Delta$ (green dotted) [119] are also shown for
comparison. Inset gives the dependence of $\Delta_{\textrm{c}}$ on
$s_{\textrm{c}}^{2}$.

We have studied above the consequences of $\widehat{P}_{\textrm{c}}$ on a few
quantities. We now investigate to which extent the limit
$\widehat{P}_{\textrm{c}}\lesssim 0.374$ holds. In order to bound
$P/\varepsilon=\widehat{P}/\widehat{\varepsilon}$ (notice that the hatted
pressure and energy density are both scaled with $\varepsilon_{\textrm{c}}$),
we need to take three generalizations of $\widehat{P}_{\textrm{c}}\lesssim
0.374$ obtained from Eq. (1) by studying the following questions:

  1. (a)

How does the $\widehat{P}/\widehat{\varepsilon}$ behave at a finite
$\widehat{r}$ for the maximum-mass configuration $M_{\textrm{NS}}^{\max}$?

  2. (b)

How does the limit $\widehat{P}_{\textrm{c}}\lesssim 0.374$ modify when we
consider stable NSs on the M-R curve away from the maximum-mass configuration?

  3. (c)

By combining (a) and (b), how does the $P/\varepsilon$ behave for stable NSs
at finite distances $\widehat{r}$ away from their centers?

For the first question, since the pressure $\widehat{P}$ and
$\widehat{\varepsilon}$ are decreasing functions of $\widehat{r}$, i.e.,

$\displaystyle\widehat{P}\approx$
$\displaystyle\widehat{P}_{\textrm{c}}+b_{2}\widehat{r}^{2}<\widehat{P}_{\textrm{c}},$
(45) $\displaystyle\widehat{\varepsilon}\approx$
$\displaystyle\widehat{\varepsilon}_{\textrm{c}}+a_{2}\widehat{r}^{2}=1+s_{\textrm{c}}^{-2}b_{2}\widehat{r}^{2}<1=\widehat{\varepsilon}_{\textrm{c}},$
(46)

we obtain by taking their ratio

$\displaystyle{P}/{\varepsilon}={\widehat{P}}/{\widehat{\varepsilon}}\approx$
$\displaystyle{\widehat{P}_{\textrm{c}}}/{\widehat{\varepsilon}_{\textrm{c}}}+\left(1-\frac{\widehat{P}_{\textrm{c}}}{s_{\textrm{c}}^{2}}\right)b_{2}\widehat{r}^{2}\approx\widehat{P}_{\textrm{c}}+\frac{1}{4}b_{2}\widehat{r}^{2}<\widehat{P}_{\textrm{c}}.$
(47)

Generally, $1-\widehat{P}_{\textrm{c}}/s_{\textrm{c}}^{2}>0$ and the
approximation $s_{\textrm{c}}^{2}\approx 4\widehat{P}_{\textrm{c}}/3$ is used
for small $\widehat{P}_{\textrm{c}}$ in the last step. This means not only
$\widehat{P}$ and $\widehat{\varepsilon}$ decrease for finite $\widehat{r}$,
but also does their ratio $\widehat{P}/\widehat{\varepsilon}$. Thus for NSs at
the maximum-mass configuration of the M-R curves, we have
$\widehat{P}/\widehat{\varepsilon}\leq\widehat{P}_{\textrm{c}}\lesssim 0.374$.
Considering the second question and for stable NSs on the M-R curve, one has
$\Psi\sim\mathrm{d}M_{\textrm{NS}}/\mathrm{d}\varepsilon_{\textrm{c}}>0$ and
Eq. (1) should be modified to

$s_{\textrm{c}}^{2}=\widehat{P}_{\textrm{c}}\left(1+\frac{1+\Psi}{3}\frac{1+3\widehat{P}_{\textrm{c}}^{2}+4\widehat{P}_{\textrm{c}}}{1-3\widehat{P}_{\textrm{c}}^{2}}\right).$
(48)

Consequently, the upper limit for $\widehat{P}_{\textrm{c}}$ constrained by
$s_{\textrm{c}}^{2}\leq 1$ should be smaller than $0.374$ due to the
positiveness of $\Psi$. Furthermore, for the last question (c), the inequality
(47) still holds and is modified for small $\widehat{P}_{\textrm{c}}$,

$\widehat{P}/\widehat{\varepsilon}\approx\widehat{P}_{\textrm{c}}+\frac{1+\Psi}{4+\Psi}b_{2}\widehat{r}^{2}<\widehat{P}_{\textrm{c}}.$
(49)

Combining the above three aspects, we find that
$P/\varepsilon=\widehat{P}/\widehat{\varepsilon}\leq\widehat{P}_{\textrm{c}}\lesssim
0.374$ holds for all stable NSs along the M-R curve either at or being close
to their centers. Nevertheless, the validity of this conclusion is still
limited to small $\widehat{r}$ (like the SSS) due to the perturbtive nature of
the expansions of $\widehat{P}(\widehat{r})$ and
$\widehat{\varepsilon}(\widehat{r})$. Whether $P/\varepsilon$ could exceed
such upper limit at even larger distances away from the centers depends on the
joint analysis of $s^{2}$ and $P/\varepsilon$, by including more higher order
contributions of the expansions.

The upper bound $P/\varepsilon\lesssim 0.374$ (at least near the NS centers)
is an intrinsic property of the TOV equations, which embody the strong-field
aspects of gravity (GR). In this sense, there is no guarantee a prior that
this bound is consistent with all microscopic nuclear EOSs. This is mainly
because the latter were conventionally constructed without considering the
strong-field ingredients of gravity, see, e.g., the sizable tension between
our causality boundary on the M-R diagram of FIG. 7 (solid black line) and
predictions of a few nuclear EOSs. The robustness of such upper bound for
$P/\varepsilon$ can be checked only by observable astrophysical
quantities/processes involving strong-field aspects of gravity such as NS M-R
data, NS-NS mergers and/or NS-black hole mergers [197].

We use the perturbative expansion of $P/\varepsilon$ and Eq. (41) to explain
the existence of a peak in the derivative part $-t$ of the SSS [119]. Dividing
Eq. (VII.2) by $\mu=\widehat{\varepsilon}-1$ gives the ratio $P/\varepsilon$
to order $\mu^{2}$ as
$P/\varepsilon\approx\widehat{P}_{\textrm{c}}-t_{\textrm{c}}\mu+(D+t_{\textrm{c}})\mu^{2}$
where $t_{\textrm{c}}=\widehat{P}_{\textrm{c}}-s_{\textrm{c}}^{2}$, this
expansion generalizes Eq. (47). Consequently, the derivative part $-t$ of
$s^{2}$ is obtained as $-t(\widehat{\varepsilon})=-t(\mu)\approx-
t_{\textrm{c}}+(2D+t_{\textrm{c}})\mu+(3J-D-t_{\textrm{c}})\mu^{2}$ (to order
$\mu^{2}$), from which the peak could be found,

$\mu_{\textrm{pk}}^{(-t)}=\widehat{\varepsilon}_{\textrm{pk}}^{(-t)}-1=\frac{1}{2}\frac{2D+t_{\textrm{c}}}{D+t_{\textrm{c}}-3J}.$
(50)

Taking for example $\widehat{P}_{\textrm{c}}\approx 0.16\mbox{$\sim$}0.30$,
this peak is shown to be located at about $-\mu_{\textrm{pk}}^{(-t)}\approx
10\%\mbox{$\sim$}5\%$ (see discussions given after Eq. (VII.2) for the
estimate on $D$ and $J$). On the other hand, using the expansion
$P/\varepsilon$ just given here, we deduce that
$t_{\textrm{c}}/2(D+t_{\textrm{c}})>0$, i.e., there would be no peak in the
non-derivative part of $s^{2}$ or this part is an increasing function of
$\varepsilon$ [119]. The value of $-t$ at $\mu_{\textrm{pk}}^{(-t)}$ is,

$-t_{\textrm{pk}}\equiv-t(\mu_{\textrm{pk}}^{(-t)})=\frac{4D^{2}+12Jt_{\textrm{c}}-3t_{\textrm{c}}^{2}}{4D-12J+4t_{\textrm{c}}}.$
(51)

Then, we have $-t_{\textrm{pk}}\approx 0.15\mbox{$\sim$}0.35$ and the non-
derivative part (of $s^{2}$), namely $P/\varepsilon\approx
0.14\mbox{$\sim$}0.28$ for $\widehat{P}_{\textrm{c}}\approx
0.16\mbox{$\sim$}0.30$. Moreover, comparing $\mu_{\textrm{pk}}^{(-t)}$ of Eq.
(50) and $\mu_{\textrm{pk}}=-D/3J$ of $s^{2}(\mu)$ (see Eq. (43)) leads us to

$\mu_{\textrm{pk}}-\mu_{\textrm{pk}}^{(-t)}=\frac{2D^{2}}{9J^{2}}\frac{1+t_{\textrm{c}}/D+3Jt_{\textrm{c}}/2D^{2}}{1-D/3J-t_{\textrm{c}}/3J}>0.$
(52)

This means the location of the peak in $s^{2}$ (i.e., $\mu_{\textrm{pk}}$)
occurs at a higher energy density than that of the peak in $-t$ (i.e.,
$\mu_{\textrm{pk}}^{(-t)}$), see also FIG. 2 of Ref. [119]. Equivalently, the
peak in $s^{2}$ is closer to the NS center than the peak in $-t$. When higher
order terms of $\mu$ are included in the expansions, the specific location of
$\mu_{\textrm{pk}}$ and $\mu_{\textrm{pk}}^{(-t)}$ may vary, however the
relation $\mu_{\textrm{pk}}-\mu_{\textrm{pk}}^{(-t)}>0$ would not be changed.

Figure 18: The $\widehat{P}_{\textrm{c}}$-dependence of coefficients
$\mathcal{C}_{\Phi}$ appeared in $\Phi\approx
1+\mathcal{C}_{\Phi}\widehat{r}^{2}+\mathcal{O}(\widehat{r}^{4})$ for
$\Phi=\widehat{P}/\widehat{P}_{\textrm{c}}$,
$\widehat{\varepsilon}/\widehat{\varepsilon}_{\textrm{c}}$,
$(\widehat{P}/\widehat{\varepsilon})/(\widehat{P}_{\textrm{c}}/\widehat{\varepsilon}_{\textrm{c}})$
and $\rho/\rho_{\textrm{c}}$.

In FIG. 18, we show the $\widehat{P}_{\textrm{c}}$-dependence of coefficient
$\mathcal{C}_{\Phi}$ appeared in the approximation $\Phi\approx
1+\mathcal{C}_{\Phi}\widehat{r}^{2}+\mathcal{O}(\widehat{r}^{4})$ for
$\Phi=\widehat{P}/\widehat{P}_{\textrm{c}},~{}\widehat{\varepsilon}/\widehat{\varepsilon}_{\textrm{c}},~{}(\widehat{P}/\widehat{\varepsilon})/(\widehat{P}_{\textrm{c}}/\widehat{\varepsilon}_{\textrm{c}})$
and $\rho/\rho_{\textrm{c}}$, for the $M_{\textrm{NS}}^{\max}$ configuration.
Here besides the relations for $\widehat{P}$, $\widehat{\varepsilon}$ and
$\widehat{P}/\widehat{\varepsilon}$, we also have the corresponding relation
for baryon density $\rho$ as (to order $\widehat{r}^{4})$,

$\widehat{\rho}\equiv\rho/\rho_{\textrm{c}}\approx
1+\left(\frac{b_{2}/s_{\textrm{c}}^{2}}{1+\widehat{P}_{\textrm{c}}}\right)\widehat{r}^{2}+\frac{1}{1+\widehat{P}_{\textrm{c}}}\left(a_{4}-\frac{{b_{2}^{2}}/{2s_{\textrm{c}}^{2}}}{1+\widehat{P}_{\textrm{c}}}\right)\widehat{r}^{4},$
(53)

obtained using the thermodynamic relation
$\rho\partial\varepsilon/\partial\rho=P+\varepsilon$ with $\rho_{\textrm{c}}$
the baryon number density at NS center. Here, $a_{4}$ characterizes the model
dependence of $\rho/\rho_{\textrm{c}}$. Without any surprise,
$\rho<\rho_{\textrm{c}}$ for finite $\widehat{r}\neq 0$. By taking
$\widehat{P}_{\textrm{c}}\approx 0.24$ for PSR J0740+6620 [115], we then
obtain numerically $\widehat{P}/\widehat{P}_{\textrm{c}}\approx
1-1.5\widehat{r}^{2}$,
$\widehat{\varepsilon}/\widehat{\varepsilon}_{\textrm{c}}\approx
1-0.8\widehat{r}^{2}$,
$(\widehat{P}/\widehat{\varepsilon})/(\widehat{P}_{\textrm{c}}/\widehat{\varepsilon}_{\textrm{c}})\approx
1-0.7\widehat{r}^{2}$, and $\rho/\rho_{\textrm{c}}\approx
1-0.6\widehat{r}^{2}$, here
$\widehat{\varepsilon}_{\textrm{c}}=\varepsilon_{\textrm{c}}/\varepsilon_{\textrm{c}}=1$.
Two reference values for $\widehat{P}_{\textrm{c}}$ at 0.24 and 0.374 are
marked using the black dash-dotted lines. It is clearly shown that all these
$\mathcal{C}_{\Phi}$’s are definitely negative for
$0\leq\widehat{P}_{\textrm{c}}\lesssim 0.374$, which are different from the
coefficient $\mathcal{C}_{\Phi}$ for $\Phi=s^{2}/s_{\textrm{c}}^{2}$. The
latter may be either positive or negative (FIG. 12 and FIG. 13). By recovering
the physical radius $r$, one can find from Eq. (45) the (reduced) pressure
$P/P_{\textrm{c}}\approx 1-0.78(r/R_{\max})^{2}$ (neglecting the intercept
0.64 in Eq. (6)), holding for perturbatively small $r/R_{\max}$. Combing (46)
and (53) gives (to order $\mu^{2}$ and $\widehat{P}_{\textrm{c}}^{2}$),

$\displaystyle\rho/\rho_{\textrm{c}}\approx$ $\displaystyle 1+\chi\mu-
s_{\textrm{c}}^{2}\chi^{2}\mu^{2}\approx\widehat{\varepsilon}-\mu\left(1+\frac{4}{3}\mu\right)\widehat{P}_{\textrm{c}}\left(1-\widehat{P}_{\textrm{c}}\right),$
(54)

where $\chi=1/(1+\widehat{P}_{\textrm{c}})$ and
$\mu=\widehat{\varepsilon}-1<0$, i.e.,
$\widehat{\rho}\approx\widehat{\varepsilon}$ to leading-order near the center
and $\widehat{\rho}\gtrsim\widehat{\varepsilon}$ considering finite
$\widehat{P}_{\textrm{c}}$. Relation (54) gives us an estimate on the location
of the possible peak in $s^{2}(\rho)$ as a function of $\rho$, e.g., by
considering $\widehat{\varepsilon}_{\textrm{pk}}-1=\mu_{\textrm{pk}}=-D/3J$ of
Eq. (43) and FIG. 14 we find that the peak is at
$\widehat{\rho}_{\textrm{pk}}=\rho_{\textrm{pk}}/\rho_{\textrm{c}}\approx
95\%$ by taking $D\approx-0.45$ and $J\approx-2.7$ (relevant for PSR
J0740+6620), i.e., it is very close to the center (the
finite-$\widehat{P}_{\textrm{c}}$ correction of (54) is $\lesssim 1\%$). The
expansion (53) is useful and we give in Appendix D an order of magnitude
estimate for the coefficient $a_{4}$ via the decreasing feature of
$\widehat{\rho}$ as a function of $\widehat{r}$.

In addition, since $P/\varepsilon$ is a decreasing function for small
$\widehat{r}$ as given in Eq. (47), the conformal anomaly
$\Delta=1/3-P/\varepsilon$ essentially increases with $\widehat{r}$ (i.e.,
when going outward) not so far ($\widehat{r}$ being small). This means that it
is disadvantage for the NS matter to become (nearly) conformal when going away
from the center not so far. This is also reflected in the behaviors of
$\gamma=s^{2}/(P/\varepsilon)$ and
$\Theta=[\Delta^{2}+(P/\varepsilon-s^{2})^{2}]^{1/2}$ at finite $\widehat{r}$
or finite $\widehat{\varepsilon}$ (or equivalently $\mu$). In fact, we can
show straightforwardly that,

$\displaystyle\gamma/\gamma_{\textrm{c}}\approx$ $\displaystyle
1+\frac{b_{2}}{s_{\textrm{c}}^{2}}\left(1+\frac{2D}{s_{\textrm{c}}^{2}}-\frac{s_{\textrm{c}}^{2}}{\widehat{P}_{\textrm{c}}}\right)\widehat{r}^{2}$
$\displaystyle\approx$ $\displaystyle
1-\frac{3D}{16\widehat{P}_{\textrm{c}}^{2}}\widehat{r}^{2},$ (55)
$\displaystyle\Theta/\Theta_{\textrm{c}}\approx$ $\displaystyle
1+\frac{b_{2}}{s_{\textrm{c}}^{2}}\frac{3t_{\textrm{c}}(1+3s_{\textrm{c}}^{2}-6\widehat{P}_{\textrm{c}}-6D)}{1+9s_{\textrm{c}}^{2}-6\widehat{P}_{\textrm{c}}(1+3s_{\textrm{c}}^{2})+18\widehat{P}_{\textrm{c}}^{2}}\widehat{r}^{2}$
$\displaystyle\approx$ $\displaystyle 1+\frac{1-6D}{8}\widehat{r}^{2},$ (56)

where
$b_{2}=-6^{-1}(1+3\widehat{P}_{\textrm{c}}^{2}+4\widehat{P}_{\textrm{c}})$ and
the coefficients in front of $\widehat{r}^{2}$ in $\gamma/\gamma_{\textrm{c}}$
and $\Theta/\Theta_{\textrm{c}}$ are both positive, the second line for each
quantity keeps only the leading-order term in $\widehat{P}_{\textrm{c}}$.
Numerically, then $\Theta/\Theta_{\textrm{c}}\approx 1+3.8\widehat{r}^{2}$ and
$\gamma/\gamma_{\textrm{c}}\approx 1+2.3\widehat{r}^{2}$ for
$\widehat{P}_{\textrm{c}}\approx 0.24$ using $D\approx-0.45$ (panel (b) of
FIG. 13). For small $\widehat{P}_{\textrm{c}}\to 0$, the second line of Eq.
(56) gives approximately $\Theta/\Theta_{\textrm{c}}\approx
1+\widehat{r}^{2}/8$ since $D\sim\widehat{P}_{\textrm{c}}^{2}\to 0$ indicated
by Eq. (55) (see also panel (b) of FIG. 13). Similarly, the $\mu$-dependence
of $\gamma$ and $\Theta$ could be obtained,

$\displaystyle\gamma/\gamma_{\textrm{c}}\approx$ $\displaystyle
1+\left(\frac{t_{\textrm{c}}}{\widehat{P}_{\textrm{c}}}+\frac{2D}{s_{\textrm{c}}^{2}}\right)\mu$
$\displaystyle\approx$ $\displaystyle
1+\frac{3D}{2\widehat{P}_{\textrm{c}}}\mu,$ (57)
$\displaystyle\Theta/\Theta_{\textrm{c}}\approx$ $\displaystyle
1+\frac{3t_{\textrm{c}}+9s_{\textrm{c}}^{2}(t_{\textrm{c}}+2D)-18\widehat{P}_{\textrm{c}}(t_{\textrm{c}}+D)}{1+9s_{\textrm{c}}^{2}-6\widehat{P}_{\textrm{c}}(1+3s_{\textrm{c}}^{2})+18\widehat{P}_{\textrm{c}}^{2}}\mu$
$\displaystyle\approx$ $\displaystyle 1+(6D-1)\widehat{P}_{\textrm{c}}\mu.$
(58)

The coefficients in front of $\mu$ in Eqs. (57) and (58) are both negative,
e.g., we have $\gamma/\gamma_{\textrm{c}}\approx 1-2.9\mu$ and
$\Theta/\Theta_{\textrm{c}}\approx 1-4.8\mu$ for
$\widehat{P}_{\textrm{c}}\approx 0.24$. Therefore, Eq. (55)$\sim$Eq. (58)
together generalize the conclusion of FIG. 16: If the matter at centers of
massive NSs was not conformal, it is likely its nearby surroundings are also
not conformal. Combining the leading-order terms of Eq. (57) and Eq. (58) or
those of Eq. (55) and Eq. (56) leads us to the correlation between the
$\gamma$ and $\Theta$ parameters as,

$\Delta\Theta\approx\left(4-{2D}/{3}\right)\widehat{P}_{\textrm{c}}^{2}\cdot\Delta\gamma,$
(59)

which holds near NS centers. Eq. (59) tells that two positive quantities
$\Delta\Theta\equiv\Theta/\Theta_{\textrm{c}}-1$ and
$\Delta\gamma\equiv\gamma/\gamma_{\textrm{c}}-1$ is positively correlated
($4-2D/3>0$ according to panel (b) of FIG. 13).

Another feature of Eq. (48) is that for NSs close to the maximum-mass point
$M_{\textrm{NS}}^{\max}$ along the M-R curve, one has
$\mathrm{d}^{2}M_{\textrm{NS}}/\mathrm{d}\varepsilon_{\textrm{c}}^{2}<0$ and
therefore the criterion (31) still holds. However, the deterministic term
$u_{\textrm{c}}=2s_{\textrm{c}}^{2}-3\widehat{P}_{\textrm{c}}$ in the
coefficient $D$ of Eq. (30) should be modified to (for small
$\widehat{P}_{\textrm{c}}$),

$u_{\textrm{c}}=2s_{\textrm{c}}^{2}-3\widehat{P}_{\textrm{c}}\approx-\frac{1-2\Psi}{3}\widehat{P}_{\textrm{c}},~{}~{}\Psi>0,$
(60)

by expanding Eq. (48) over $\widehat{P}_{\textrm{c}}$ as
$s_{\textrm{c}}^{2}\approx(4+\Psi)\widehat{P}_{\textrm{c}}/3$ (Newtonian
limit). A positive $\Psi$ tends to make the correction $D$ positive and
therefore reduce the probability of
$s_{\textrm{c}}^{2}<s^{2}(\widehat{\varepsilon})$ (see Eq. (23)). For example,
a canonical NS is less likely to have the crossover in the core than a massive
NS (assuming they have similar radii), since the former is more likely (than
the latter) to be on climbing stage of the M-R curve. Our analysis is
consistent with Ref. [152] which predicted the increasing of $s^{2}$ is
reduced even to disappear when going out as the NS mass decreases from
$M_{\textrm{NS}}^{\max}$ to $1.4M_{\odot}$. On the other hand, for massive NSs
with similar radii of PSR J0740+6620 (thus the $\widehat{P}_{\textrm{c}}$’s
are also similar), the reduction on the probability of $D<0$ due to a (small)
positive $\Psi$ is expected to be small as they are near the maximum-mass
configuration [189, 190, 191].

##  IX Distinguishing Compactness from Stiffness in Neutron Stars, and pQCD
prediction from General Relativity requirement on Trace Anomaly

With both analytical analyses and numerical examples, we have investigated
separately in detail in the previous two sections the radial variation of SSS
and the various bounds on the trace anomaly. To this end, it is useful to
summarize the main differences between the compactness and stiffness of NSs.
It is also important to emphasize the respective applicabilities of the pQCD
prediction and GR limit on trace anomaly in NSs.

Figure 19: Qualitative sketch of the reduction of $\widehat{P}$,
$\widehat{\varepsilon}$, $\widehat{P}/\widehat{\varepsilon}$ and
$\widehat{\rho}$ (these quantities characterize the compactness/denseness of
the matter) to order $\widehat{r}^{2}$ and the probable increasing of $s^{2}$
(characterizing the stiffness of the matter) near NS centers. Here
$b_{2}=-6^{-1}(1+3\widehat{P}_{\textrm{c}}^{2}+4\widehat{P}_{\textrm{c}})$ and
the coefficient $D$ depends on $\widehat{P}_{\textrm{c}}$ (see panel (b) of
FIG. 13).

According to the perturbative expressions (45)$\sim$(47) (as well as FIG. 18),
the energy density $\varepsilon$, the pressure $P$, their ratio
$P/\varepsilon$ and the baryon density $\rho$ (see Eq. (53)) are all
decreasing functions of $\widehat{r}$ definitely. While on the other hand, the
SSS $s^{2}$ has sizable probabilities to be enhanced when going outward from
centers of NSs, see sketches shown in FIG. 19. This demonstrates the
fundamental difference between the stiffness (characterized effectively by
$s^{2}$) and the denseness/compactness (characterized effectively by $\rho$,
$P$, $\varepsilon$ or $P/\varepsilon$), although they are closely related to
each other (see, e.g., FIG. 8 or the $s_{\textrm{c}}^{2}$ of Eq. (1) itself).
An illustrative/useful example is that $s_{\textrm{c}}^{2}$ (of Eq. (1))
approaches $1/3$ (stiffness) earlier than $\widehat{P}_{\textrm{c}}\to 1/3$
(compactness). Physically, it is because the EOS of NS matter (especially near
the centers) is nonlinear (otherwise $s^{2}=P/\varepsilon$, $\gamma=1$ and
$\Theta=[\Delta^{2}+(P/\varepsilon-s^{2})^{2}]^{1/2}=|\Delta|$).

Closely related is the trace anomaly (or the conformality measure)
$\Delta=1/3-P/\varepsilon$ [119, 198]. We sketch in FIG. 20 the $\Delta$ as a
function of energy density $\varepsilon$, here $\varepsilon_{0}\approx
150\,\textrm{MeV}/\textrm{fm}^{3}$ is the fiducial energy density at
$\rho_{\textrm{sat}}$, around which the low-energy nuclear theories constrain
the $\Delta$ quite well. The pQCD bound $\Delta=0$ applies effectively at very
large energy densities $\varepsilon\gtrsim 50\varepsilon_{0}\approx
7.5\,\textrm{GeV}/\textrm{fm}^{3}$ [119, 199], which is far larger than the
energy density reacheable in NSs. Thus, it is possibly relevant but not
fundamental for explaining the observed $P/\varepsilon\gtrsim 1/3$ in massive
NSs (using microscopic/phenomenological models). On the other hand, we have
demonstrated that a GR bound on
$\widehat{P}_{\textrm{c}}=P_{\textrm{c}}/\varepsilon_{\textrm{c}}$ and
$P/\varepsilon$ (near NS centers) naturally emerges when dissecting
perturbatively the TOV equations without using any specific EOS model [115].
In this sense, the GR bound on $\Delta$ (with the $\varepsilon$ being roughly
around (4$\sim$8$)\varepsilon_{0}$) is likely more relevant/fundamental than
the pQCD prediction for nuclear EOS in NSs, although the latter may influence
the extraction of $\Delta$ [193]. Whether these two specific bounds (being
very close to each other) are inner-related may deepen our understanding on
the connection between GR and the microscopic theories of elementary
particles. Our above discussions and findings are not out of bounds as
ultimately properties of compact objects are determined by the Hamilton’s
principle using the total action of the whole system including gravity, matter
(including nuclear matter, dark matter and energy) and their couplings [200].

Figure 20: Sketch of the patterns for $\Delta=1/3-P/\varepsilon$ in NSs. The
$\Delta$ is well constrained around the fiducial density
$\varepsilon_{0}\approx 150\,\textrm{MeV}/\textrm{fm}^{3}$ by low-energy
nuclear theories and is predicted to vanish due to conformality of the matter
at $\varepsilon\gtrsim 50\varepsilon_{0}$ using pQCD theories. Massive and
compact NSs provide a unique opportunity to probe the negativeness of $\Delta$
in certain energy density regions, where the GR bound on $\Delta$ (with
$\varepsilon$ being around ($4\mbox{$\sim$}8)\varepsilon_{0}$) is expected to
be more relevant than the pQCD bound ($\varepsilon\gtrsim 50\varepsilon_{0}$)
in these NSs.

As a negative $\Delta$ is unlikely to be observed in ordinary NSs (e.g., NSs
with masses $\sim 1.7M_{\odot}$), see FIG. 17, the evolution of $\Delta$ is
probably like the green line in FIG. 20. An (unconventional) exception may
come from light but very compact NSs, e.g., a $1.7M_{\odot}$ NS with radius
about 9.3 km (therefore $\varepsilon_{\textrm{c}}\approx
1.86\,\textrm{GeV}/\textrm{fm}^{3}$ and $P_{\textrm{c}}\approx
654\,\textrm{MeV}/\textrm{fm}^{3}$ and $\widehat{P}_{\textrm{c}}\approx
0.351$) has its $\Delta_{\textrm{c}}\approx-0.02$ (FIG. 17). On the other
hand, massive and compact NSs (masses $\gtrsim 2M_{\odot}$) have the most
relevance to observe a negative $\Delta$ (as indicated by the magenta line in
FIG. 20), and how the negative $\Delta$ evolves to the pQCD bound may tell
more on the properties of supradense matter. Unfortunately, the region with
$\varepsilon\gtrsim 8\varepsilon_{0}$ is largely inaccessible in NSs (sketched
by the light-cyan rectangle band) due to their self-gravitating nature (see
FIG. 9 which indicates that observations of massive NSs should inevitably put
an upper bound on the central energy density as
$\varepsilon_{\textrm{c}}\lesssim\varepsilon_{\textrm{ult}}$).

##  X Summary

Using the central SSS $s_{\textrm{c}}^{2}$ and NS mass/radius scaling obtained
from analyzing perturbatively structures of the scaled TOV equations [115], we
studied the radial variation of SSS, trace anomaly and several closely related
properties of NSs in an EOS-model independent manner. We found that the
reduced pressure
$\widehat{P}_{\textrm{c}}=P_{\textrm{c}}/\varepsilon_{\textrm{c}}$ is the most
relevant quantity determining the onset of a continuous crossover of SSS in
the cores of massive NSs. With sizable probabilities the crossover occurs near
the centers of massive NSs if $\widehat{P}_{\textrm{c}}\lesssim 0.3$, e.g.,
with a probability $\gtrsim 63\%$ for $\widehat{P}_{\textrm{c}}\approx 0.24$
(PSR J0740+6620). The resulting peak of $s^{2}$ or its derivative part
(defined as $s^{2}-P/\varepsilon$) near NS centers is generally demonstrated.

With the help of the universal correlations of
$M_{\textrm{NS}}^{\max}$-$\Gamma_{\textrm{c}}$ and
$R_{\max}$-$\nu_{\textrm{c}}$ [115] and the nonlinear dependence of
$s_{\textrm{c}}^{2}$ on $\widehat{P}_{\textrm{c}}$ (of Eq. (1)), a new
causality boundary for NS M-R curve as $R_{\max}/\textrm{km}\gtrsim
4.73M_{\textrm{NS}}^{\max}/M_{\odot}+1.14$ was obtained. It is shown to be
excellently consistent with several NS mass/radius observations and puts a
more stringent constraint on the dense matter EOS compared to the ones
available in the literature. Moreover, the NS maximum compactness parameter is
limited to $\lesssim 0.313\cdot(1-1.14\,\textrm{km}/R_{\max})\lesssim 0.283$
(using $R_{\max}\approx 12\,\textrm{km}$). The difference between $s^{2}\leq
1$ (stiffness) and $P/\varepsilon\leq 1$ (compactness) in NSs is clarified,
with its physical origin being traced back to the nonlinear characteristic of
the core EOS in NSs. While $s^{2}\leq 1$ is more relevant than
$P/\varepsilon\leq 1$ in NSs, the latter is shown to be $\lesssim 0.374$ or
equivalently $\Delta=1/3-P/\varepsilon\gtrsim-0.041$ around the centers of
stable NSs. Such bound is a direct consequence of GR (encapsulated in the TOV
equations), and is expected to be more fundamental than the pQCD prediction of
dense matter EOS for NSs. Closely related, the matter in cores of massive NSs
is unlikely to be conformal if the criteria
$\Theta=[\Delta^{2}+(P/\varepsilon-s^{2})^{2}]^{1/2}\lesssim 0.2$ and
$\gamma=\mathrm{d}\ln P/\mathrm{d}\ln\varepsilon\lesssim 1.75$ were adopted.
Moreover, the violation of conformal bound is found to basically depend on
$\widehat{P}_{\textrm{c}}$, e.g., the bound is broken if
$M_{\textrm{NS}}^{\max}\gtrsim 1.9M_{\odot}$ as $s_{\textrm{c}}^{2}$ increases
with $M_{\textrm{NS}}^{\max}/M_{\odot}$ (under the assumption NSs with masses
about $1.3\mbox{$\sim$}2.3M_{\odot}$ have similar radii $\approx$12 km).
Furthermore, the ultimate energy density $\varepsilon_{\textrm{ult}}$ and
pressure $P_{\textrm{ult}}$ allowed in NSs are estimated. The observations of
massive NSs may induce corresponding upper limits for
$\varepsilon_{\textrm{ult}}$ and $P_{\textrm{ult}}$, e.g., the existence of a
2.08$M_{\odot}$ NS leads to $\varepsilon_{\textrm{ult}}\lesssim
1.32\,\textrm{GeV}/\textrm{fm}^{3}$ and $P_{\textrm{ult}}\lesssim
494\,\textrm{MeV}/\textrm{fm}^{3}$, respectively. This improves previous
constraints on the same quantities.

Finally, unlike most existing studies of NSs in the literature, our analyses
are carried out without using any model EOS for NS matter thus largely free of
the uncertainties in modeling the dense matter EOS. We have demonstrated that
some properties of the underlying NS EOS (pressure versus energy density) can
be inferred directly from the observational data alone. This in turn can help
constrain predictions of dense matter EOS based on nuclear many-body theories.
Overall, our work contributes to realizing the ultimate goal of understanding
properties of strong-field gravity, dense matter and their couplings in
compact objects.

## Acknowledgements

We would like to thank Peng-Cheng Chu for helpful discussions. This work is
supported in part by the U.S. Department of Energy, Office of Science, under
Award No. DE-SC0013702, the CUSTIPEN (China-U.S. Theory Institute for Physics
with Exotic Nuclei) under US Department of Energy Grant No. DE-SC0009971, the
National Natural Science Foundation of China under Grant No.12235010.

##  Appendix A  Proof of Eq. (30) for a General $K$

From Eq. (28), we can solve for $d_{3}$ in terms of other $d_{k}$’s (using the
property that $s^{2}=0$ for $\widehat{\varepsilon}=0$),

$d_{3}=-2\widehat{P}_{\textrm{c}}+s_{\textrm{c}}^{2}-\sum_{k=4}^{K}(k-2)d_{k},$
(A61)

where we split the summation $\sum_{k=3}^{K}(k-2)d_{k}$ in Eq. (28) into
$d_{3}+\sum_{k=4}^{K}(k-2)d_{k}$. Therefore, according to the basic definition
of $D$ of Eq. (22), we obtain,

$\displaystyle D=$
$\displaystyle\sum_{k=1}^{K}\frac{k(k-1)}{2}d_{k}=d_{2}+3d_{3}+\sum_{k=4}^{K}\frac{k(k-1)}{2}d_{k}$
$\displaystyle=$ $\displaystyle
2s_{\textrm{c}}^{2}-3\widehat{P}_{\textrm{c}}+\sum_{k=4}^{K}\left(\frac{k(k-1)}{2}-(k-1)-(k-2)\right)d_{k}$
$\displaystyle=$ $\displaystyle
2s_{\textrm{c}}^{2}-3\widehat{P}_{\textrm{c}}+\sum_{k=4}^{K}\frac{(k-2)(k-3)}{2}d_{k},$
(A62)

here expression for $d_{2}$ of Eq. (29) is used. In summary, condition
$d_{1}=0$ and the sum rules of (21) together enables us to rewrite $D$ in the
form of Eq. (30), and for the special case of $K=3$ the $D$ becomes
deterministic.

##  Appendix B Explaining the Sign of $D$ and $J$

In the main text, we analytically study the case of $K=3$ (see FIG. 11), where
a sharp change for the probability of $D<0$ to $D>0$ occurs at
$\widehat{P}_{\textrm{c}}\approx 0.105$. For more general values of $K$,
similar phenomenon may happen. In this appendix, we try to explain why the
coefficient $D$ tends to be negative (positive) for small (large)
$\widehat{P}_{\textrm{c}}$, using the expansions with $K=4$. We have two
relations, i.e., $d_{2}+d_{3}+d_{4}=\widehat{P}_{\textrm{c}}$ and
$2d_{2}+3d_{3}+4d_{4}=s_{\textrm{c}}^{2}$, from which we can solve for $d_{2}$
and $d_{3}$ to be expressed in terms of $d_{4}$. The general condition $0\leq
s^{2}\leq 1$ now reads $0\leq
2d_{2}\widehat{\varepsilon}+3d_{3}\widehat{\varepsilon}^{2}+4d_{4}\widehat{\varepsilon}^{3}\leq
1$. Putting the expressions of $d_{2}$ and $d_{3}$ into it gives the upper and
lower limit for $d_{4}$, i.e., $d_{4}^{\textrm{(l)}}\leq d_{4}\leq
d_{4}^{\textrm{(u)}}$, where

$\displaystyle d_{4}^{\textrm{(l)}}=$
$\displaystyle\frac{\widehat{\varepsilon}(\widehat{P}_{\textrm{c}}-s_{\textrm{c}}^{2})-\widehat{P}_{\textrm{c}}+\widehat{\varepsilon}^{-1}}{4\widehat{\varepsilon}^{2}-3\widehat{\varepsilon}-1},~{}~{}d_{4}^{\textrm{(u)}}=\frac{\widehat{\varepsilon}(\widehat{P}_{\textrm{c}}-s_{\textrm{c}}^{2})-\widehat{P}_{\textrm{c}}}{4\widehat{\varepsilon}^{2}-3\widehat{\varepsilon}-1}.$
(B63)

Adding the deterministic term
$u_{\textrm{c}}=-3\widehat{P}_{\textrm{c}}+2s_{\textrm{c}}^{2}$ of Eq. (32),
we obtain correspondingly
$D^{\textrm{(u/l)}}=u_{\textrm{c}}+d_{4}^{\textrm{(u/l)}}$ (see Eq. (30)).

Figure b1: Explanation on why $D$ tends to be negative (positive) for small
(large) $\widehat{P}_{\textrm{c}}$, using $K=4$. See text for details.

Figure b2: The same as FIG. b1 but for the coefficient $J$ of Eq. (38).

The results are summarized in FIG. b1. Here in each panel, the magenta (blue)
line is for $D^{(\textrm{u})}$ ($D^{\textrm{(l)}}$), and the orange line marks
the general boundary of $2D<\sigma_{\textrm{c}}^{2}s_{\textrm{c}}^{2}$ (see
the estimate algorithm of (31)). The positive and negative regions for $D$ are
shown by the grey and lavender bands, respectively. Moreover, $D\leq
D^{(\textrm{u})}$ indicates that $D$ should be smaller than the lowest point
$D_{1}$ on the magenta line (solid green diamond in panel (c)). Similarly,
$D\geq D^{\textrm{(l)}}$ implies that $D$ should be larger than the highest
point $D_{2}$ on the blue line (shallow red circle in panel (c)). When
$\widehat{P}_{\textrm{c}}$ is small (such as panels (a) and (b)), the
magnitude of red circle is larger than that of the greed diamond, indicating
that $D$ is more probable to be negative than be positive. As
$\widehat{P}_{\textrm{c}}$ increases one eventually encounters $D_{1}\approx-
D_{2}$ (as panel (c) shows), $D_{2}\approx 0$ and $D_{2}\gtrsim 0$ (panel
(d)). The situation $D_{1}\approx-D_{2}$ roughly marks the value of
$\widehat{P}_{\textrm{c}}$ corresponding to equal probability of $D>0$ and
$D<0$, this is about $\widehat{P}_{\textrm{c}}\approx 0.3$ (as verified by
simulation given in the main text, see FIG. 12). When
$\widehat{P}_{\textrm{c}}$ increases even further, the probability of $D<0$
(of $D>0$) quickly decreases (increases). The coefficient $J$ could be
similarly analyzed, i.e.,
$J=3d_{4}+(s_{\textrm{c}}^{2}-\widehat{P}_{\textrm{c}})/3$ where
$d_{4}^{\textrm{(l)}}\leq d_{4}\leq d_{4}^{\textrm{(u)}}$ (see (38) and
(B63)). Therefore $J^{\textrm{(l)}}\leq J\leq J^{\textrm{(u)}}$, and we
probably have $J<0$ though the allowed region of $J<0$ (of $J>0$) eventually
shrinks (expands) as $\widehat{P}_{\textrm{c}}$ increases, see FIG. b2. For
other values of $K$, the general features are similar as FIG. b1 and FIG. b2
and are consistent with the results given in the main text (see FIG. 12 and
FIG. 13).

Figure b3: The same as FIG. b1 but adopting the Newtonian limit
$s_{\textrm{c}}^{2}\approx 4\widehat{P}_{\textrm{c}}/3$ (left panel) and
conformal limit $s_{\textrm{c}}^{2}\approx\widehat{P}_{\textrm{c}}$ (right
panel). The principle of causality requires then
$\widehat{P}_{\textrm{c}}\lesssim 3/4$ and $\widehat{P}_{\textrm{c}}\lesssim
1$ for these two situations, respectively. Symbols and lines have the same
meaning of FIG. b1.

Moreover, if one artificially adopts the Newton approximation for
$s_{\textrm{c}}^{2}$ as $s_{\textrm{c}}^{2}\approx
4\widehat{P}_{\textrm{c}}/3$ (panel (a) of FIG. b3) or the conformal limit
$\gamma_{\textrm{c}}=s_{\textrm{c}}^{2}/\widehat{P}_{\textrm{c}}\approx
1\leftrightarrow s_{\textrm{c}}^{2}\approx\widehat{P}_{\textrm{c}}$ (panel (b)
of FIG. b3), then even a large $\widehat{P}_{\textrm{c}}\approx 0.6$ may still
tend to induce a negative $D$ (compared with FIG. b1). In particular, the $D$
in conformal limit is very likely negative. In addition, if
$s_{\textrm{c}}^{2}$ is nearly a constant (independent of
$\widehat{P}_{\textrm{c}}$), then the condition
$D<\sigma_{\textrm{c}}^{2}s_{\textrm{c}}^{2}/2$ alone tells that $D<0$ since
$\sigma_{\textrm{c}}^{2}=\mathrm{d}s_{\textrm{c}}^{2}/\mathrm{d}\widehat{P}_{\textrm{c}}=0$
(see the inequality (26)). In this case if $s_{\textrm{c}}^{2}=1$, then its
nearby $s^{2}(\widehat{\varepsilon})$ keeps 1, otherwise
$s^{2}(\widehat{\varepsilon})>s_{\textrm{c}}^{2}$. These results clearly show
that the nonlinear dependence of $s_{\textrm{c}}^{2}$ on
$\widehat{P}_{\textrm{c}}$ (of Eq. (1)) is fundamental to account for the
eventual change on the sign of $D$, and therefore it may influence the
possible continuous crossover occurring in cores of NSs.

##  Appendix C Perturbative Corrections to $s_{\textmd{c}}^{2}$

In this appendix, we estimate the first a few perturbative corrections to
$s_{\textrm{c}}^{2}$ from the $b_{4}$-term. The scheme is given as follows: We
first notice that [115],

$b_{4}=-\frac{1}{2}b_{2}\left(\widehat{P}_{\textrm{c}}+\frac{4+9\widehat{P}_{\textrm{c}}}{15s_{\textrm{c}}^{2}}\right),$
(C64)

where
$b_{2}=-6^{-1}(1+3\widehat{P}_{\textrm{c}}^{2}+4\widehat{P}_{\textrm{c}})$