###### Abstract

Recent inference results of the sound velocity in the cores of neutron stars
are summarized. Implications for the equation of state and the phase structure
of highly compressed baryonic matter are discussed. In view of the strong
constraints imposed by the heaviest known pulsars, the equation of state must
be very stiff in order to ensure the stability of these extreme objects. This
required stiffness limits the possible appearance of phase transitions in
neutron star cores. In view of these empirical findings, much discussed issues
such as the quest for a phase transition towards restored chiral symmetry and
the active degrees of freedom in cold and dense baryonic matter, are explored.

###### keywords:

dense baryonic matter; neutron stars; phase transitions

1 1 0 https://doi.org/ Constraints on Phase Transitions in Neutron Star Matter
Title Len Brandes* and Wolfram Weise* Correspondence<EMAIL_ADDRESS><EMAIL_ADDRESS>
##  1 Introduction

Within recent decades the observational data base for neutron stars has
expanded substantially. Its detailed analysis sharpened the empirical
constraints on the equation of state (EoS) of matter deep inside the cores of
these extreme objects. In particular, the existence of heavy neutron stars
with masses around and above two solar masses implies that the EoS, i.e.
pressure $P(\varepsilon)$ as function of energy density $\varepsilon$, must be
sufficiently stiff in order to support such massive compact stars against
gravitational collapse. Some previously discussed simple forms of exotic
matter could therefore be excluded when their corresponding EoS’s turned out
to be too soft.

The quest for a possible phase transition in cold and dense baryonic matter,
from hadronic to quark degrees of freedom, has been a topic of prime interest
for a long time. In this context a summary will be given of the present state
of knowledge about the neutron star EoS as inferred directly from the existing
observational data. Of particular interest is the speed of sound in the cores
of neutron stars. Its behaviour as a function of energy density is a sensitive
indicator for phase transitions or continuous changes (crossovers) of degrees
of freedom in the star’s composition.

Primary sources of information for the inference of the EoS are neutron star
masses deduced from Shapiro delay measurements Demorest2010 ; Antoniadis2013 ;
Fonseca2016 ; Arzoumanian2018 ; Cromartie2020 ; Fonseca2021 , combined
determinations of masses and radii inferred from X-ray data detected by the
NICER telescope Riley2019 ; Miller2019 ; Riley2021 ; Miller2021 ; Salmi2022 ,
and gravitational wave signals from binary neutron star mergers observed by
the LIGO and Virgo Collaborations Abbott2018 ; Abbott2019 ; Abbott2020 . These
data have served as basic input for a large variety of Bayesian inference
analyses Raaijmakers2021 ; Pang2021 ; Legred2021 ; Biswas2022 ; Ecker2022 ;
Altiparmak2022 ; Huth2022 ; Lim2022 ; Annala2023 ; Somasundaram2023 ;
Essick2023 ; Brandes2023 ; Mroczek2023 ; Brandes2023a ; Han2023  in search for
a range of data-compatible equations of state with controlled uncertainties.
Low-density constraints from nuclear physics around the equilibrium density of
nuclear matter, $\rho_{0}=0.16\,$fm-3, are commonly introduced by referring to
results of chiral effective field theory (ChEFT) Drischler2022 . At
asymptotically high densities, perturbative QCD (pQCD) is applicable and
provides a constraint to be matched in extrapolations far beyond the
conditions realized in neutron star cores Gorda2021 ; Komoltsev2022 . The
Bayesian inference studies have recently been further extended Brandes2023a ,
given the supplementary observational data reported in Romani2022  (the black
widow pulsar PSR J0952-0607 with a mass $M=2.35\pm 0.17\,M_{\odot}$ in solar
mass units).

This report presents a survey of results obtained in the detailed inference
analysis of Ref.Brandes2023a  and discusses possible interpretations. A
central topic is the evaluation of constraints and possible evidence for a
first-order phase transition within the posterior uncertainties of the EoS
deduced from the existing empirical neutron star data. Later parts of this
manuscript refer to ongoing considerations about the structure and composition
of dense matter in neutron star cores. In view of the constraints emerging
from the analysis of the observational data, especially of the heaviest known
neutron stars, we discuss the quest for a chiral phase transition, quark-
hadron continuity and its realisation as a soft crossover, a Fermi liquid
approach to neutron star matter, and other related issues.

##  2 Neutron star matter equation of state

###  2.1 Observational constraints

####  2.1.1 Neutron star masses and radii

Information about the masses of the heaviest neutron stars derives primarily
from precise Shapiro time delay measurements of pulsars orbiting in binary
systems with white dwarf companions. Three such massive objects, PSR
J1614–2230 Demorest2010 ; Fonseca2016 ; Arzoumanian2018 , PSR J0348+0432
Antoniadis2013  and PSR J0740+6620 Cromartie2020 ; Fonseca2021 , have been
established in the past:

PSR J1614–2230 $\displaystyle M=1.908\pm 0.016\,M_{\odot}~{},~{}$ (1) PSR
J0348+0432 $\displaystyle M=2.01\pm 0.04\,M_{\odot}~{},$ (2) PSR J0740+6620
$\displaystyle M=2.08\pm 0.07\,M_{\odot}~{}.$ (3)

Recently the heaviest neutron star observed so far has been reported
Romani2022 , a black widow pulsar,

PSR J0952-0607 $\displaystyle M=2.35\pm 0.17\,M_{\odot}~{}.~{}$ (4)

This result was obtained based on the detection of the black widow’s companion
star. PSR J0952-0607 is also one of the fastest rotating pulsars. With a spin
period of 1.4 ms it requires radius-dependent corrections for rotational
effects as described in Brandes2023a . For example, the rotating mass of a
$2.35\,M_{\odot}$ star with an assumed radius of $R=12$ km decreases by 3% to
an equivalent non-rotating mass of about 2.28 $M_{\odot}$.

Radii of neutron stars together with their masses can be inferred from X-ray
profiles of rotating hot spot patterns measured with the NICER telescope. Two
neutron stars have been investigated in this way Riley2019 ; Riley2021 :

PSR J0030+0451 $\displaystyle\qquad
M=1.34^{+0.15}_{-0.16}\,M_{\odot}~{},\qquad$ $\displaystyle
R=12.71^{+1.14}_{-1.19}\,\text{km}~{},$ (5) PSR J0740+6620
$\displaystyle\qquad M=2.072^{+0.067}_{-0.066}\,M_{\odot}~{},\qquad$
$\displaystyle R=12.39^{+1.30}_{-0.98}\,\text{km}~{}.$ (6)

All of the mass and radius data are included as input in setting constraints
for the neutron star EoS. In the case of PSR J0740+6620 the analysis of the
NICER measurement includes the Shapiro delay result. Only the data from NICER
have therefore been used in the Bayesian inference data set. An alternative
analysis of the NICER data has also been performed by a second independent
group Miller2019 ; Miller2021 . Both results are mutually compatible within
their uncertainties.

####  2.1.2 Binary neutron star mergers and tidal deformabillities

Gravitational wave signals produced by the merger of two neutron stars in a
binary have been detected by the LIGO and Virgo collaborations Abbott2019 ;
Abbott2018 ; Abbott2020 . These signals are interpreted using theoretical
waveform models that depend on the mass ratio of the two neutron stars,
$M_{2}/M_{1}$, and a mass-weighted combination of their tidal deformabilities.
Based on information from the GW170817 event the (dimensionless) tidal
deformability for a 1.4 solar mass neutron star has been deduced Abbott2018 :

$\displaystyle\Lambda_{1.4}=190^{+390}_{-120}~{}.$ (7)

The GW170817 event was further investigated together with electromagnetic
signals. In one such analysis Fasano2019  the following masses and tidal
deformabilities were reported for the individual neutron stars in the binary :

$\displaystyle\qquad M_{1}=1.46^{+0.13}_{-0.09}\,M_{\odot}~{},\qquad$
$\displaystyle\Lambda_{1}=255^{+416}_{-171}~{},$ (8) $\displaystyle\qquad
M_{2}=1.26^{+0.09}_{-0.12}\,M_{\odot}~{},\qquad$
$\displaystyle\Lambda_{2}=661^{+858}_{-375}~{}.$

The event GW190425 which could have also been a binary neutron star merger,
was also included in the inference analysis Abbott2020 .

###  2.2 Inference of sound velocity and EoS in neutron stars

Based on the data in the previous subsection the squared sound velocity and
the equation of state,

$\displaystyle c_{s}^{2}(\varepsilon)=\frac{\partial
P(\varepsilon)}{\partial\varepsilon}~{}~{}~{}\text{and}~{}~{}~{}P(\varepsilon)=\int_{0}^{\varepsilon}d\varepsilon^{\prime}\,c_{s}^{2}(\varepsilon^{\prime})~{},$
(9)

of neutron star matter can be inferred using Bayesian methods. In our recent
work Brandes2023 ; Brandes2023a  the ChEFT constraint at low baryon density is
used in the Bayes inference procedure as a likelihood (not as a prior) within
a conservative window of baryon densities Essick2020 , $\rho\lesssim
1.3\,\rho_{0}$. The equations of state are interpolated to pQCD results at
asymptotically high densities Komoltsev2022 .

Figure 1: Squared speed of sound Brandes2023a  as function of energy density
inferred from the empirical data set listed in Section 2.1. Median (solid
curve) and posterior credible bands at levels of 68% (grey band) and 95%
(dashed curves) are displayed. Intervals are indicated for 68% ranges of
central energy densities in the cores of 1.4 and 2.3 $M_{\odot}$ neutron
stars.

The result of the inferred sound velocity in Figure 1 shows, within the $68\%$
credible bands, a rapid increase of $c_{s}^{2}$ beyond the conformal limit of
1/3 at energy densities relevant for a wide range of neutron stars with masses
$M\sim$ 1.4 - 2.3 $M_{\odot}$. A quantification in terms of a corresponding
Bayes factor Legred2021 ; Brandes2023 ; Brandes2023a  demonstrates indeed
extreme evidence for $c_{s}^{2}$ exceeding the conformal bound inside neutron
stars.

Figure 2 presents the inferred posterior bands for the EoS. Notably, the
additional inclusion of the heaviest pulsar PSR J0952-0607 in the data set
leads to a more rapidly rising pressure as compared to previous inference
results which do not incorporate this data. Within the range of energy
densities 0.5 - 0.8 GeV/fm3 reached in the core of a 2.3 solar mass star and
corresponding to baryon densities between 3.0$\rho_{0}$ and 4.6$\rho_{0}$, the
pressure even exceeds that of the time-honored APR EoS APR1998 .

Figure 2: Equation of state $P(\varepsilon)$ Brandes2023a  deduced from the
inferred sound velocity based on the empirical data set listed in Section 2.1.
Median, 68% and 95% posterior credible bands are displayed as in Figure 1. The
APR EoS APR1998  (dotted line) is shown for comparison. Also shown for
orientation is the baryon density scale, $\rho/\rho_{0}$ (in units of the
equilibrium density of nuclear matter $\rho_{0}=0.16$ fm-3), computed for the
median of $P(\varepsilon)$.

A further quantity of interest is the chemical potential associated with the
conserved baryon number,

$\displaystyle\mu(\rho)=\frac{\partial\varepsilon(\rho)}{\partial\rho}=\frac{P+\varepsilon}{\rho}~{},$
(10)

where the latter equality follows from the Gibbs-Duhem relation at zero
temperature, $P+\varepsilon=\mu\rho$. This baryon chemical potential can be
computed for any given equation of state within the uncertainty bands in
Figure 2. The resulting $\mu(\rho)$ bands in Figure 3 show a steep rise at
baryon densities $\rho>2\,\rho_{0}$, indicating strongly repulsive
correlations at work. The agnostic approach for the inference of the sound
speed on which these results are based does not permit to distinguish between
different species of constituents inside neutron stars, so $\mu$ represents
the total chemical potential from all active degrees of freedom carrying
baryon number: $\mu=\sum_{i}x_{i}\mu_{i}$, where $x_{i}=\rho_{i}/\rho$ is the
fraction corresponding to each baryonic species $i$. But irrespective of
whether these species are nucleons, other baryonic composites or quarks, the
empirically inferred behaviour of $\mu$ demonstrates that these degrees of
freedom must be correlated by strongly repulsive forces as the density
increases, in order to build up the necessary high pressure in the cores of
neutron stars to support two solar masses and beyond. We return to this
discussion again at a later stage.

Figure 3: Baryon chemical potential $\mu(\rho)$ in neutron star matter,
normalized at $\mu(\rho=0)=939$ MeV. Shown are posterior credible bands
Brandes2023a  at 68% level (grey band) and 95% level (dashed lines). 2.

The approach to conformality in strongly interacting matter at high baryon
density has been an important theme in several recent investigations
Fujimoto2022 ; Marczenko2023 ; Rho2023 ; Ma2023 . A key quantity is the trace
anomaly measure,

$\displaystyle\Delta=\frac{g_{\mu\nu}T^{\mu\nu}}{3\varepsilon}=\frac{1}{3}-\frac{P(\varepsilon)}{\varepsilon}~{}.$
(11)

Conformal matter is characterized by an equation of state $P=\varepsilon/3$,
so the signature for its occurance is $\Delta\rightarrow 0$ over an extended
energy density range. The posterior credible bands for the trace anomaly
measure deduced from the data-based EoS are shown in Figure 4. Starting from
$\Delta=1/3$ at zero energy density, the median of $\Delta$ turns negative at
$\varepsilon\sim 0.7\,\text{GeV/fm}^{3}$, entering a high-pressure domain with
$P>\varepsilon/3$. This crossing appears within the range of energy densities
possibly realized in the cores of the heaviest neutron stars. At much higher
energy densities beyond those displayed in Figure 3, pQCD then implies a
switch back to positive $\Delta$ before the asymptotic limit
$\Delta\rightarrow 0$ is approached.

Figure 4: Median and credible bands Brandes2023a  of the trace anomaly
measure, $\Delta=\frac{1}{3}-\frac{P}{\varepsilon}$, at 68% level (grey band)
and 95% level (dashed lines).

###  2.3 Selected neutron star properties

Numerically solving the Tolman-Oppenheimer-Volkov (TOV) equations leads to
posterior credible bands for the mass-radius relation of neutron stars as
displayed in Figure 5. Notably, the median of $M(R)$ suggests an almost
constant radius $R$, independent of $M$.

Figure 5: Posterior credible bands Brandes2023a  of the mass-radius relation
at 68% level (grey band) and 95% level (dashed line), compared to the NICER
data analysis by Riley et al. of PSR J0030+0451 and PSR J0740+6620 Riley2019 ;
Riley2021 . In addition the 68% mass-radius credible interval of the
thermonuclear burster 4U 1702-429 (blue) is displayed Nattila2017  (which was
not included in the Bayesian analysis).

The density profiles $\rho(r)$ computed for neutron stars in the mass range
1.4 - 2.3 $M_{\odot}$ reach central baryon densities, $\rho_{c}=\rho(r=0)$,
that are systematically below $5\,\rho_{0}$ in the cores of even the heaviest
stars. For example,

$\displaystyle\rho_{c}(1.4\,M_{\odot})=(2.6\pm
0.4)\,\rho_{0}\qquad\text{and}\qquad\rho_{c}(2.3\,M_{\odot})=(3.8\pm
0.8)\,\rho_{0}~{},$ (12)

at the 68% level Brandes2023a . In a baryonic picture of neutron stars, this
implies that the average distance between baryons even in the highly
compressed star center exceeds 1 fm, more than twice the characteristic short-
range hard-core distance of 1/2 fm in the nucleon-nucleon interaction. We
return to this point at a later stage.

Figure 6: Posterior credible bands Brandes2023a  of the tidal deformability as
a function of neutron star mass, at 68% level (grey band) and 95% level
(dashed lines), compared to values Abbott2019 ; Fasano2019  deduced from the
merger event GW170817.

For each given EoS, the TOV equations can be solved in combination with
equations for the tidal deformability $\Lambda$. This leads to posterior
probability bands, $\Lambda(M)$, as a function of neutron star mass as shown
in Figure 6. The large uncertainties of the empirical values deduced from the
analysis of the gravitational wave events do so far not allow to constrain the
tidal deformability bands to a stronger degree.

##  3 Constraints on phase transitions in neutron stars

The nature and location of the transition from dense nuclear or neutron matter
to quark and gluon degrees of freedom is still largely unknown. Lattice
computations have been successful in exploring the QCD phase structure at high
temperatures and vanishing baryon chemical potential, establishing that the
transition is a continuous crossover at a temperature $T_{c}\simeq 155$ MeV.
However, at nonzero baryon densities the notorious sign problem in the
Euclidean action prevents a similarly systematic approach. As a consequence,
studies of the phase diagram at low temperatures and high baryon densities are
primarily based on models with varying assumptions about the relevant active
degrees of freedom Sumiyoshi2023 .

A key element of this discussion is the question about a possible phase
transition from spontaneously broken to restored chiral symmetry in QCD. At
high baryon densities beyond $\rho>2\,\rho_{0}$ various hypotheses have been
explored, ranging from a chiral first-order phase transition to continuous
hadron-quark crossovers Fukushima2011 ; Holt2016 ; McLerran2019 ; Baym2019 ;
Fukushima2020 ; Kojo2022 . The behaviour of the sound speed $c_{s}$ is a prime
indicator for phase transitions or crossovers. This is schematically
illustrated in Figure 7 which shows typical patterns of such transitions as
they would show up in the sound speed. For example, a first-order phase
transition with Maxwell construction is characterized by a region of constant
pressure over an interval of density (or energy density) in which two phases
coexist. Between the lower and upper end points of this phase coexistence
interval the squared sound velocity $c_{s}^{2}=\partial P/\partial\varepsilon$
jumps to zero and back. In a first-order phase transition with Gibbs
construction the pressure in the mixed phase is not constant but $\partial
P/\partial\varepsilon$ still changes discontinuously Han2019 . A crossover has
no such discontinuities but features a pronounced maximum in $c_{s}^{2}$.
Purely nucleonic scenarios lead instead to a continuously rising sound speed.

Figure 7: Characteristic behaviours of the squared sound velocity in the
presence of a phase transition or a crossover.

###  3.1 Evidence against a very low squared sound speed in neutron stars

One indication for the possible appearance of a first-order phase transition
would be a very low squared sound velocity $(c_{s}^{2}\leq 0.1)$ within the
range of energy densities relevant for neutron stars. Let us introduce
$c_{s,min}$ as the minimum speed of sound following a maximum at lower
densities. Given the inferred bands of probability distributions for
$c_{s}^{2}$ in Section 2.2 and Figure 1, one can quantify the evidence for two
competing scenarios: minimum $c^{2}_{s,min}>0.1$ versus $c^{2}_{s,min}\leq
0.1$, in terms of corresponding Bayes factors, ${\cal
B}^{c^{2}_{s,min}>0.1}_{c^{2}_{s,min}\leq 0.1}$. A detailed analysis of such
Bayes factors has been carried out in Brandes2023a  for the minimum speed of
sound up to different neutron star maximum masses, $M_{max}$. Bayes factors
exceeding 10 or 100 indicate strong or extreme evidence, respectively, for the
scenario where the minimum $c_{s}^{2}$ remains larger than 0.1. This is the
case for all masses $M_{max}\leq 2.1\,M_{\odot}$ as demonstrated in Table 3.1.
Only for even heavier neutron stars does this evidence turn out to be lower.
It is worth noting that the inclusion of the heavy black widow pulsar PSR
J0952-0607 plays an important role strengthening the evidence against
$c_{s}^{2}\leq 0.1$ further in comparison with previous analyses which do not
inlude this object Brandes2023 . Squared sound speeds $c_{s}^{2}>0.1$ in the
cores of neutron stars with $M_{max}\leq 2.0\,M_{\odot}$ were also found in
Altiparmak2022 ; Annala2023  without inclusion of the black widow pulsar.

[H] Bayes factors ${\cal B}^{c^{2}_{s,min}>0.1}_{c^{2}_{s,min}\leq 0.1}$
comparing EoS samples with competing scenarios: a) minimum squared speed of
sound (following a maximum) with $c^{2}_{s,min}>0.1$, versus b) equations of
state with $c^{2}_{s,min}\leq 0.1$. The minimum speeds of sound are computed
up to maximum neutron star masses as indicated. Taken from Brandes2023a .
$M_{max}/M_{\odot}$ 1.9 2.0 2.1 2.2 2.3 ${\cal
B}^{c^{2}_{s,min}>0.1}_{c^{2}_{s,min}\leq 0.1}$ 500.9 229.8 15.0 3.6 2.2 These
limiting bounds on very low speeds of sound suggest that the appearance of a
strong first-order phase transition with Maxwell construction is unlikely for
most of the known neutron stars. An additional feature of a Maxwell-
constructed first-order transition is an extended phase coexistence region.
This will be examined more closely in the following subsection.

###  3.2 Evidence against a strong first-order phase transition in the cores
of neutron stars

The relatively moderate baryon densities inferred in the cores of neutron
stars together with the evidence against very small sound speeds can be
complemented by an additional study that is specific to first-order phase
transitions with Maxwell construction. Such phase transitions are
characterized by a domain of phase coexistence that extends over a certain
range of baryon densities, $\Delta\rho$. The width of this domain,
$\Delta\rho/\rho$ (with $\rho$ the density at which the coexistence interval
starts), is a measure of the strength of the phase transition, i.e. the
magnitude of surface tension between the two coexisting phases. We refer to a
‘strong’ first-order phase transition if the width of the mixed phase is
$\Delta\rho/\rho>1$. One example is the liquid-gas phase transition in
symmetric nuclear matter at temperatures $T<15$ MeV Wellenhofer2014 ;
Brandes2021 . In contrast a ‘weak’ first-order phase transition has
$\Delta\rho/\rho$ small compared to unity.

Figure 8: Illustration of the constraint on the maximum width of a Maxwell-
constructed coexistence region for a first-order phase transition within the
68% credible band of $P(\varepsilon)$ (see Figure 2) following ref.
Brandes2023a .

The posterior credible bands inferred from neutron star data, as displayed in
Figure 2, permit a systematic study of the maximum possible phase coexistence
widths over the range of relevant energy densities. Starting from
$P(\varepsilon)$ the Gibbs-Duhem relation is used to reexpress pressure as a
function of baryon density, $P(\rho)$. In this way posterior credible bands
for $P(\rho)$ can be derived. These constrain the maximum possible widths of
mixed phase domains, $(\Delta\rho/\rho)_{\text{max}}$. An example is shown in
Figure 8. In the analysis Brandes2023a  it turns out that these maximum
possible phase coexistence regions are narrow:
$(\Delta\rho/\rho)_{\text{max}}\simeq 0.2$ at the 68% level. Even within the
95% credible bands this width does not exceed 0.3. In fact one finds that
$(\Delta\rho/\rho)_{\text{max}}$ stays nearly constant as a function of baryon
density $\rho$ (taken at the starting point of the mixed-phase interval) over
the whole region of densities $\rho\simeq 2-5\,\rho_{0}$ relevant for neutron
stars.

###  3.3 Intermediate summary

The conclusion drawn in the reported analysis is therefore that only a weak
first-order phase transition with a mixed-phase width $\Delta\rho/\rho\leq
0.2-0.3$ in a Maxwell construction can be realized inside neutron stars within
the posterior credible bands of the inferred equation of state. A strong
first-order phase transition is thus very unlikely to occur Brandes2023a . It
is worth pointing out that this conclusion is entirely based on the analysis
of observational data, independent of any specific EoS model. As will be
discussed later, this rules out some relativistic mean-field models which
suggested the occurance of a strong (chiral) first-order phase transition
already at densities below those encounterd in neutron star cores. On the
other hand the empirical uncertainty bands still leave room for a continuous
transition such as a hadrons-to-quarks crossover Baym2019 ; McLerran2019 ;
Fukushima2020 ; Kojo2022 . The possiblity of a purely baryonic equation of
state with monotonically rising speed of sound is also not excluded
Brandes2023 .

##  4 Phenomenology and models

The previous sections have shown how present and forthcoming neutron star data
(masses, radii, tidal deformabilities) provide constraints for the equation of
state and for possible phase transitions in cold, dense baryonic matter. The
enhanced stiffness of the EoS at high baryon densities is a necessary feature:
sufficiently high pressures are required to support neutron stars with masses
around and beyond $2\,M_{\odot}$ at relatively large radii ($R\simeq 11-13$
km). The central core densities are not as extreme as was previously imagined.
Even in the core of a 2.3 $M_{\odot}$ neutron star the baryon density at its
center stays below five times the density of equilibrium nuclear matter (at
68% credible level, see Eq. (12)). The upper border of the 95% credible
interval of the inferred central density
$\rho_{c}(2.3\,M_{\odot})/\rho_{0}=3.8^{+1.6}_{-1.3}$ Brandes2023a  lies
slightly above $5\,\rho_{0}$, but with very low probability. This observation
has consequences for the interpretation of the possible structure and
composition of neutron star matter that we discuss in the following.

###  4.1 Reminder of low-energy nucleon structure and a two-scales scenario

Spontaneouly broken chiral symmetry as the long-wavelength manifestation of
QCD governs the low-energy structure and dynamics of hadrons including
nucleons and pions. As chiral Nambu-Goldstone bosons, pions play a
distinguished role in this context. Models based on chiral symmetry often view
the nucleon as a complex system of two scales TW2001 : a compact hard core
that hosts the three valence quarks and thus incloses the baryon number, and a
surrounding quark-antiquark cloud in which pions figure prominently as the
‘soft’ degrees of freedom.

Such a two-scales scenario is manifest in empirical form factors and sizes of
the nucleon. Consider for example the proton and neutron electromagnetic form
factors and their slopes at zero momentum transfer which determine the
corresponding mean squared radii. The empirical r.m.s. proton charge radius,
$\langle r_{p}^{2}\rangle^{1/2}=0.840\pm 0.003$ fm, combined with the slope of
the neutron electric form factor, $\langle r_{n}^{2}\rangle=-0.105\pm 0.006$
fm2, gives the isoscalar and isovector mean-square radii of the nucleon,
$\langle r^{2}_{S,V}\rangle=\langle r_{p}^{2}\rangle\pm\langle
r_{n}^{2}\rangle$ with resulting values Lin2022 :

$\displaystyle\sqrt{\langle r_{S}^{2}\rangle}\simeq
0.78\,\text{fm}~{},\quad~{}~{}\sqrt{\langle r_{V}^{2}\rangle}\simeq
0.90\,\text{fm}~{}.$ (13)

Each of the nucleon form factors related to a current with index $i$ has a
representation in terms of a dispersion relation,

$\displaystyle
G_{i}(q^{2})=G_{i}(0)+\frac{q^{2}}{\pi}\int_{t_{0}}^{\infty}dt\,\frac{\text{Im}\,G_{i}(t)}{t(t-q^{2}-i\epsilon)}~{}.$
(14)

Mean-square radii are given as:

$\displaystyle\langle
r_{i}^{2}\rangle=6\frac{dG_{i}(q^{2})}{dq^{2}}\Big{|}_{q^{2}=0}=\frac{6}{\pi}\int_{t_{0}}^{\infty}\frac{dt}{t^{2}}S_{i}(t)~{},$
(15)

where the distribution $S_{i}(t)=\text{Im}G_{i}(t)$ represents the spectrum of
intermediate hadronic states through which the external probing field couples
to the respective nucleon current. For example, the isovector charge radius
reflects the interacting two-pion cloud of the nucleon governed by the $\rho$
meson and a low-mass tail extending down to the $\pi\pi$ threshold,
$t_{0}=4m_{\pi}^{2}$. The isoscalar charge radius is related to the three-pion
spectrum that is strongly dominated by the narrow $\omega$ meson Kaiser2019
and starts at $t_{0}=9m_{\pi}^{2}$. The isoscalar charge form factor of the
nucleon, $G_{S}^{E}(q^{2})$ (with $G_{S}^{E}(0)=1$), is particularly suitable
to discuss a delineation between ‘core’ and ‘cloud’ parts of the nucleon
Brown1986 ; Meissner1987 . The vector meson dominance principle implies, in
its simplest version, a representation of the form:

$\displaystyle
G_{S}^{E}(q^{2})=\frac{F_{B}(q^{2})}{1+|q^{2}|/m_{\omega}^{2}}~{}.$ (16)

The form factor $F_{B}(q^{2})$ of the baryon number distribution in the
nucleon core acts as a source for the $\omega$ field that propagates with its
mass $m_{\omega}$. Introducing the mean-square radius of the baryon core,
$\langle r_{B}^{2}\rangle=6\frac{dF_{B}(q^{2})}{dq^{2}}|_{q^{2}=0}$, the mean-
square isoscalar charge radius becomes

$\displaystyle\langle r_{S}^{2}\rangle=\langle
r_{B}^{2}\rangle+\frac{6}{m_{\omega}^{2}}~{}.$ (17)

Using $m_{\omega}=783$ MeV the estimated baryonic core radius is:

$\displaystyle\sqrt{\langle r_{B}^{2}\rangle}\simeq 0.47\,\text{fm}~{}.$ (18)

A nucleon core size of about 1/2 fm is characteristic of chiral ‘core + cloud’
models. It also keeps up in more detailed treatments of the spectral
distributions governing the nucleon form factors Lin2022 . Inclusion of
additional $\phi$ meson and $\pi\rho$ continuum contributions in the spectral
function of $G_{S}^{E}(q^{2})$ move the core radius to just slightly larger
values.

Consider as another example the form factor associated with the axial vector
current of the nucleon. The corresponding mean-square axial radius deduced
from neutrino-deuteron scattering data is reported as Hill2018 :

$\displaystyle\langle r_{A}^{2}\rangle=(0.46\pm 0.22)\,\text{fm}^{2}~{}.$ (19)

A schematic axial vector dominance picture would assign a dominant part of the
‘cloud’ contribution to this form factor through the spectrum of the $a_{1}$
meson with its large width. If an approximate scale of this ‘cloud’ part is
identified with $\delta\langle r_{A}^{2}\rangle\sim 6/m_{a}^{2}$ using the
physical $a_{1}$ mass, $m_{a}\simeq 1.23$ GeV, one finds

$\displaystyle\sqrt{\langle r_{A}^{2}\rangle}_{\text{core}}=[\langle
r_{A}^{2}\rangle-\delta\langle r_{A}^{2}\rangle]^{1/2}\simeq 0.55\,\text{fm}$
(20)

with an estimated uncertainty of about 25%.

Yet another interesting piece of information is the mass radius of the proton
deduced from $J/\Psi$ photoproduction data Kharzeev2021 . It involves the form
factor of the trace, $T_{\mu}^{\mu}$, of the nucleon’s energy-momentum tensor
and is supposed to be dominated by gluon dynamics at the center of the
nucleon. Low-mass (e.g. two-pion) components are suppressed. The result quoted
in Kharzeev2021 ,

$\displaystyle\langle r^{2}\rangle_{\text{mass}}^{1/2}=(0.55\pm
0.03)\,\text{fm}~{},$ (21)

is once again remarkably close to an assumed ‘core’ size scale of $\sim 1/2$
fm.

These empirical considerations motivate a picture of the nucleon as containing
a compact ‘hard’ core in which the valence quarks and their baryon number are
confined, and a ‘soft’ surface of quark-antiquark pairs forming a mesonic
cloud. This structure has implications for the behaviour of nucleons in dense
baryonic matter. With a core size $R_{\text{core}}\sim 1/2$ fm and a cloud
range given e.g. by the proton charge radius, $R_{\text{cloud}}\sim 0.84$ fm,
there is a significant separation of volume scales in vacuum:
$(R_{\text{cloud}}/R_{\text{core}})^{3}\sim 5$.

This scale separation is expected to increase further in dense baryonic
matter, for the following reasons. The properties of the soft multi-pion cloud
are closely tied to spontaneously broken chiral symmetry and the approximate
Nambu-Goldstone boson nature of the pion. The size of this cloud is expected
to increase with density, along with the decreasing in-medium pion decay
constant, $f_{\pi}^{*}(\rho)$, which acts as a chiral order parameter. The
baryonic core, on the other hand, is governed by gluon dynamics, without a
leading connection to chiral symmetry in QCD. This core is therefore expected
to be quite stable against changes from increasing density up until the
compact hard cores begin to touch and overlap.

Figure 9: Sketch of low- and high-density baryonic matter. Baryons (e.g.
nucleons) are viewed as valence quark cores surrounded by clouds of quark-
antiquark pairs (e.g. chiral meson clouds). At densities $\rho\gtrsim
2-3\,\rho_{0}$ percolation of quark-antiquark pairs over larger distances
starts as indicated. Valence quark cores begin to touch and overlap at baryon
densities $\rho\gtrsim 5\,\rho_{0}$.

What arises in this way is a two-scales scenario for dense baryonic matter as
described in Fukushima2020  and sketched in Figure 9. At $\rho\simeq\rho_{0}$
the tails of the meson clouds of nucleon pairs overlap, resulting in two-body
exchange forces. As the average distance between nucleons decreases with
increasing density, around $\rho\gtrsim 2-3\,\rho_{0}$, the soft clouds of
$q\bar{q}$ pairs start to be delocalized. Their mobility expands over larger
distances in a percolation process involving larger numbers of nucleons. In
conventional nuclear physics terms, this corresponds to the emergence of many-
body forces, the strength of which grows with increasing density. At that
stage the baryonic cores are still separated but subject to increasingly
repulsive Pauli principle effects. The cores begin to touch and overlap at
average nucleon-nucleon distances $d\lesssim 1$ fm, corresponding to densities
$\rho\gtrsim 5\,\rho_{0}$. Further compression of baryonic matter would still
have to overcome the strongly repulsive hard core in the nucleon-nucleon
interaction. Recalling the inferred credible intervals of the central
densities in heavy neutron stars, Eq. (12), one concludes that a phase
transition to valence quark deconfinement does not seem likely in a two-scales
picture and under the conditions in the cores of neutron star, unless they are
extremely heavy.

An interesting and closely related result emerges from a detailed analysis of
$y-$scaling in electron-nucleus scattering at large momentum transfers
($|\vec{q}|\gtrsim 1$ GeV) and low energy transfers Benhar2023 . The
persistently observed $y-$scaling under these kinematical conditions implies
that the electrons still scatter from strongly correlated pairs of nucleons,
rather than quarks, at short distances. The conclusion drawn in Benhar2023  is
that at local densities as large as five times $\rho_{0}$, nuclear matter
still appears to behave as a collection of nucleons.

###  4.2 Quark-hadron continuity and crossover

While a strong first-order phase transition in neutron star matter seems
unlikely based on current empirical observations, a continuous crossover from
hadrons to quarks is still possible within the present data-driven
constraints. Such a scenario is realised for example in the QHC21 equation of
state Kojo2022 . It features a smooth interpolation between low and high
densities from nuclear to quark matter regimes. The quark matter phase is
described by a three-flavor Nambu–Jona-Lasinio (NJL) model including pairing
(diquark) degrees of freedom and a strongly repulsive vector coupling between
quarks. In order to build up the necessary pressure at high densities, this
vector coupling must be comparable to or larger than the strength of the
standard scalar-pseudoscalar interaction in the NJL model. The density at
which the interpolated turnover to quark matter takes place in the QHC21 EoS
is chosen as $\rho\gtrsim 3.5\,\rho_{0}$, within the range of central
densities that can be reached in $M\gtrsim 2\,M_{\odot}$ neutron stars. A
crossover from hadrons to quarks could also involve an intermediate phase of
quarkyonic matter McLerran2019 ; Fujimoto2023 .

###  4.3 Chiral symmetry restoration in dense matter: from first-order phase
transition to crossover

The quest for chiral symmetry restoration at high baryon densities, a
transition from the spontaneously broken Nambu-Goldstone realisation to the
unbroken Wigner-Weyl mode, has been a persistent issue for a long time. A
possible first-order chiral phase transition and the existence of a
corresponding critical end point in the QCD phase diagram have always been
topics of prime interest Stephanov1998 ; Fukushima2003 ; Fischer2014 . Early
hypotheses concerning the occurance of a first-order phase transition were
frequently based on Nambu–Jona-Lasinio type (NJL) models in mean-field
approximation Asakawa1989 ; Klimt1990 ; Scavenius2001 , later extended by
incorporating confinement aspects through added Polyakov loop degrees of
freedom Roessner2007 ; Roessner2008 ; Hell2010 .

The empirical constraints on strong first-order phase transitions in dense
neutron star matter as described in Section 3.2 are strikingly at variance
with previous mean-field (e.g. NJL model) predictions. These suggested that a
first-order chiral phase transition should already appear at relatively low
baryon densities, $\rho\sim 2-3\,\rho_{0}$. A possible explanation for this
discrepancy can be found in Drews2017 ; Brandes2021  where a chiral nucleon-
meson (ChNM) field theory was used to explore the effects of fluctuations
beyond mean-field (MF) approximation. The starting point is a relativistic
chiral Lagrangian, ${\cal L}(\Psi;\pi,\sigma;v_{\mu})$, shaped around a linear
sigma model with nucleons $(\Psi)$, pions $(\pi)$ and a scalar $(\sigma)$
field. Short-range dynamics is parametrized in terms of heavy isoscalar and
isovector vector fields $(v_{\mu})$. The expectation value
$\langle\sigma\rangle$ of the scalar field acts as a chiral order parameter
normalized in the vacuum to the pion decay constant, $f_{\pi}\simeq 92$ MeV.
Two classes of fluctuations beyond MF were then systematically studied: first,
vacuum fluctuations which introduce an additional term proportional to
$\sigma^{4}\ln(\sigma/f_{\pi})$ in the MF partition function; and secondly,
fluctuations involving pion loops and nucleon particle-hole excitations. The
vacuum fluctuations can be included in an extended mean-field (EMF)
approximation Skokov2010 . Fluctuations involving pion and nucleon loops are
computed using non-perturbative functional renormalization group (FRG)
methods. The parameters of the ChNM model, in particular those related to
short-distance dynamics, are fixed to reproduce empirical nuclear physics data
Drews2017 ; Brandes2021 .

Figure 10: Chiral order parameters in symmetric nuclear matter and neutron
matter at temperature $T=0$ as function of baryon density in units of nuclear
ground state equilibrium density, $\rho_{0}=0.16$ fm-3. Dotted lines: liquid-
gas phase transition (L-G) in symmetric nuclear matter. Dashed lines: first-
order chiral phase transitions emerging from the mean-field (MF) approximation
of a relativistic chiral nucleon-meson field-theoretical model. Solid lines
show results of extended mean-field calculations (EMF, with inclusion of
fermionic vacuum fluctuations), and functional renormalization group
calculations (FRG), based on the same ChNM model. Figures adapted from Refs.
Brandes2021 ; Drews2017 .

Figure 10 demonstrates the important role of fluctuations beyond mean-field
approximation for the chiral order parameter $\langle\sigma\rangle$. In
symmetric nuclear matter the mean-field approximation of the ChNM model
correctly reproduces the first-order liquid-gas phase transition at low
density. But the chiral order parameter also displays a strong first-order
chiral phase transition with a Maxwell-constructed phase coexistence region
starting already below $2\,\rho_{0}$. For neutron matter, which has no liquid-
gas phase transition, the MF approximation nevertheless predicts a first-order
chiral phase transition at densities around $\rho\sim 3\,\rho_{0}$, well
within the range of densities realised in neutron stars. However, in both
nuclear and neutron matter, the inclusion of fermionic vacuum fluctuations
(i.e. the effect of the ground state zero-point energy) in the extended mean-
field (EMF) approximation converts the first-order chiral phase transition
into a smooth crossover and shifts it to densities $\rho>5\,\rho_{0}$. This
effect is further enhanced by the additional fluctuations included in the full
FRG calculation as demonstrated in Figure 10. As a result, chiral symmetry
restoration is relegated to very high baryon densities beyond the inferred
core densities in even the heaviest neutron stars (see Eq. (12)). A comparable
impact of fluctuations on the phase structure is seen in alternative chiral
models Gupta2012 ; Zacci2018 .

Another approach based on chiral symmetry is the parity-doublet model. In this
model the active coupled baryonic degrees of freedom are the nucleon with
spin-parity $1/2^{+}$ and its chiral partner with spin-parity $1/2^{-}$, where
the latter is identified with the $N^{*}(1535)$ resonance. Spontaneous chiral
symmetry breaking in vacuum manifests itself in the mass splitting of these
two states, while the in-medium restoration pattern of this symmetry is
signalled by the N and N* masses becoming degenerate. A recent detailed
analysis Eser2023  of the chiral order parameter in this model using extended
mean-field approximation finds a chiral phase transition in nuclear matter,
but at extremely high densities $(\rho>10\,\rho_{0})$, far beyond the density
scales reached in neutron stars.

###  4.4 Dense baryonic matter: a Fermi liquid picture

The constraints on the equation of state inferred from the current empirical
data still permit an interpretation of neutron star core matter in terms of
baryonic degrees of freedom, such as a system dominated by neutrons APR1998
with small fractions of protons and perhaps hyperons Lonardoni2015 ;
Gerstung2020 ; Leong2023 . The inferred baryon chemical potential
$\mu(\rho)=\frac{\partial\varepsilon(\rho)}{\partial\rho}$ shown in Figure 3
does not distinguish between different species of baryons. Its behaviour
displays nonetheless increasingly strong correlations at high densities. It is
instructive to analyze the gross properties of this state of matter using the
Landau theory of relativistic Fermi liquids Baym1976 . Here we perform a
schematic study assuming "neutron-like" quasiparticles Friman2019 , with
baryon number $B=1$ and a density-dependent mass $m(\rho)$, ignoring other
small admixtures in the composition of the dense medium. These quasiparticles
are characterized by their (relativistic) Landau effective mass $m_{L}^{*}$ at
Fermi momentum $p_{F}=(3\pi^{2}\,\rho)^{1/3}$,

$\displaystyle m_{L}^{*}(\rho)=\sqrt{p_{F}^{2}+m^{2}(\rho)}~{},$ (22)

together with an effective potential, $U(\rho)$, so that the baryon chemical
potential can be written Friman2019 :

$\displaystyle\mu(\rho)=m_{L}^{*}(\rho)+U(\rho)~{}.$ (23)

The median of the $\mu(\rho)$ posterior credible bands in Figure 3 is now
taken as a guiding starting point to extract the baryonic quasiparticle
properties. With an educated ansatz for $m(\rho)$, the density dependence of
the potential $U(\rho)$ can then be deduced and further discussed. One
possible choice is to take $m(\rho)$ (with $m(0)=939$ MeV) from the non-
perturbative FRG calculation employing the chiral nucleon-meson field-
theoretical model Drews2017  discussed in Section 4.1. The resulting Landau
effective mass $m_{L}^{*}(\rho)$ is shown in Figure 11 together with the
potential $U(\rho)=\mu(\rho)-m_{L}^{*}(\rho)$. The 95% credible band of
$\mu(\rho)$ in Figure 3 leads to an uncertainty of about 15% for $U$ at high
densities. It is instructive to fit the resulting quasiparticle potential by a
series in powers of baryon density for $\rho\lesssim 5\,\rho_{0}$ (with
$\rho_{0}=0.16\,\text{fm}^{-3}$ as before):

$\displaystyle
U(\rho)=\sum_{n}u_{n}\left(\frac{\rho}{\rho_{0}}\right)^{n}~{}.$ (24)

The coefficients fitted to the median of $\mu(\rho)$ are:

$u_{1}=90.9\,\text{MeV},\quad~{}u_{2}=15.3\,\text{MeV},\quad~{}u_{3}=3.2\,\text{MeV},\quad~{}u_{4}=-0.4\,\text{MeV}~{}.$
(25)

This pattern reflects a hierarchy of many-body correlations, recalling that
the term linear in density represents two-body interactions, the term of order
$\rho^{2}$ corresponds to short-range three-body forces, and so forth. The
role of the repulsive $N$-body terms with $N>2$ is quite significant: at
$\rho\simeq 4\,\rho_{0}$, corresponding to an average distance of about 1 fm
between the baryonic quasiparticles, these terms contribute as much as the
two-body forces to the potential $U$ and generate the strong pressure to
support heavy neutron stars. Of course these statements rely on the ansatz for
the density dependent mass $m(\rho)$ which is guided by FRG calculations of
pure neutron matter. A small fraction of protons in beta-equilibrated matter
will not change this picture substantially. However, neutron star core
compositions that deviate qualitatively from this simplified picture may lead
to different conclusions.

Figure 11: Landau effective mass
$m_{L}^{*}(\rho)=\sqrt{p_{F}^{2}+m^{2}(\rho)}$ and potential $U(\rho)$ of
quasiparticles representing the median of the posterior distribution of the
baryon chemical potential $\mu(\rho)$ (see Figure 3).

Finally, consider the dimensionless Fermi liquid parameters, $F_{0}$ and
$F_{1}$, of the spin-averaged quasiparticle interaction. In terms of the
quasiparticle mass and potential, they are given as Friman2019 :

$F_{0}(\rho)=\frac{p_{F}}{\pi^{2}}\left[m(\rho)\frac{\partial
m}{\partial\rho}+m_{L}^{*}(\rho)\frac{\partial
U}{\partial\rho}\right]~{},\quad\quad~{}F_{1}(\rho)=-\frac{3U(\rho)}{\mu(\rho)}~{}.$
(26)

Further useful relations are: $1+F_{0}=N(0)(\partial\mu/\partial\rho)$ with
the density of quasiparticle states at the Fermi surface,
$N(0)=m_{L}^{*}p_{F}/\pi^{2}$, and $1+F_{1}/3=1-U/\mu=m_{L}^{*}/\mu$.

The result in Figure 12 shows a strongly increasing $F_{0}$ at high densities.
This reflects once again the growing importance of many-body correlations as
matter gets more and more compact. Such repulsive correlations are responsible
for the increase of the sound velocity beyond its canonical conformal limit,
$c_{s}>1/\sqrt{3}$, as seen in Figure 1. The Fermi liquid parameter $F_{1}$ is
smaller in magnitude and has a negative slope indicating the decreasing
effective mass at higher densities.

Figure 12: Landau Fermi liquid parameters derived from the quasiparticle
properties that are based on the median of the data-inferred baryon chemical
potential $\mu$ (see Fig.3).

The Landau parameter $F_{0}$ displays the typical behavior of a strongly
correlated Fermi system. It is interesting though to observe that in
comparison with the leading Landau parameters in liquid 3He, the correlations
in neutron star core matter are not extraordinarily strong. Values of
$F_{0}\approx 9.3$ are reported for 3He at zero pressure, and $F_{0}\approx
68.2$ at a pressure of 27 bar Baym1991 . Accordingly, the Fermi liquid
parameters for matter in the density range, $\rho\lesssim 5\,\rho_{0}$, inside
even the heaviest neutron stars, are significantly smaller. In this
perspective the dense baryonic medium encountered in the center of neutron
stars is perhaps not as extreme as often imagined.

##  5 Concluding remarks

The neutron star data base, recently extended by adding the heaviest so far
observed pulsar (PSR J0952-0607), requires a much stiffer neutron star matter
equation of state compared to previous analyses. As a consequence, baryon
densities reached in the cores of even very heavy neutron stars do not exceed
about five times the density of normal nuclear matter. In chiral model
descriptions of the nucleon which suggest a scale separation between a valence
quark core and a mesonic cloud, the compact cores with radii of about 1/2 fm
would begin to touch and overlap only in the inner cores of extremely heavy
neutron stars. The appearance of a strong first-order phase transition becomes
unlikely under such conditions. This is consistent with the empirical results
which suggest that extended phase coexistence regions and small minimum sound
speeds are not favoured by the current data. This interpretation is further
endorsed by the fact that fluctuations beyond mean-field approximation convert
a possible first-order chiral phase transition into a smooth crossover at high
densities. At the same time, a hadrons-to-quarks crossover scenario, or
alternatively, baryonic matter as a strongly correlated relativistic Fermi
liquid, remain both possible within the data-driven inferred constraints.

Research partially supported by Deutsche Forschungsgemeinschaft (DFG grant
TRR110) and National Natural Science Foundation of China (NSFC grant no.
11621131001) through the Sino-German CRC110 "Symmetries and the Emergence of
Structure in QCD", and by the DFG Excellence Cluster ORIGINS.

Not applicable.

###### Acknowledgements.

Stimulating discussions and communications with Jean-Paul Blaizot, Kenji
Fukushima, Norbert Kaiser and Mannque Rho are gratefully acknowledged. One of
the authors (L.B.) thanks the Japan Society for the Promotion of Science and
The University of Tokyo for support and hospitality during his summer 2023
International Fellowship for Research in Japan. Not applicable. References

## References

  * (1) Demorest, P.B.; Pennucci, T.; Ransome, S.M.; Roberts, M.S.E. A two-solar mass neutron star measured using Shapiro delay. Nature 2010, 467, 1081–1083.
  * (2) Antoniadis, J.; Freire, P.C.C.; Wex, N.; Tauris, T.M.; Lynch, R.S.; van Kerkwijk, M.H.; Kramer, M.; Bassa, C.; Dhillon, V.S.; Driebe, T.; et al.. A massive pulsar in a compact relativistic binary. Science 2013, 340, 1233232, 448–502.
  * (3) Fonseca, E.; Pennucci, T.; Ellis, J.A. ; Stairs, I.H.; Nice, D.J.; Ransom, S.M.; Demorest, P.B.; Arzoumanian, Z.; Crowther, K. ; Dolch, T.; et al. The NANOGrav nine-year data set: mass and geometric measurements of binary millisecond pulsars. Astrophys. J. 2016, 832, 167–189.
  * (4) Arzoumanian, Z.; Brazier, A; Burke-Spolaor, S; Chamberlin, S.; Chatterjee, S.; Christy, B.; Cordes, J.M.; Cornish, N.J.; Crawford, F.; Cromartie, H.T.; et al. The NANOGrav 11-year data set: high-precision timing of 45 millisecond pulsars. Astrophys. J. Suppl. 2018, 235, 37–78.
  * (5) Cromartie, H.T.; Fonseca, E.; Ransom, S.M.; Demorest, P.B.; Arzoumanian, Z.; Blumer, P.R.; Brook, P.R.; DeCesar, M.E.; Dolch, T.; Ellis, J.A.; et al. Relativsitic Shapiro delay measurements of an extremely massive millisecond pulsar. Nature Astron. 2020, 4, 72–76 .
  * (6) Fonseca, E.; Cromartie, H.T.; Pennucci, T.T.; Ray, P.S.; Kirichenko, A.Y.; Ransom, S.M.; Demorest, P.B.; Stairs, I.H.; Arzoumanian, Z.; Guillemot, L.; et al. Refined mass and geometric measurements of the high-mass PSR J0740+6620. Astrophys. J. Lett. 2021 915, L12, 15pp.
  * (7) Riley, T.E.; Watts, A.L.; Bogdanov, S.; Ray, P.S.; Ludlam, R.M.; Guillot, S.; Arzoumanian, Z.; Baker, C.L.; Bilous, A.V.; Chakrabarty, D.; et al. A NICER view of PSR J0030+0451: millisecond pulsar parameter estimation. Astrophys. J. Lett. 2019, 887, L21, 60pp.
  * (8) Miller, M.C.; Lamb, F.K.; Dittmann, A.J.; Bogdanov, S.; Arzumanian, Z.; Gendreau, K.C.; Guillot, S.; Harding, A.K.; Ho, W.C.G.; Lattimer, J.M. ; et al. PSR J0030+0451 mass and radius from NICER data and implications for the properties of neutron stars. Astrophys. J. Lett. 2019, 887, L24, 28pp.
  * (9) Riley, T.E.; Watts, A.L.; Ray, P.S.; Bogdanov, S.; Guillot, S.; Morsink, S.M.; Bilous, A.V.; Arzoumanian, Z.; Choudury, D.; Daneva, J.S.; et al. A NICER view of the massive pulsar PSR J0740+6620 informed by radio timing and XMM-Newton spectroscopy. Astrophys. J. Lett. 2021, 918, L27, 30pp.
  * (10) Miller, M.C.; Lamb, F.K.; Dittmann, A.J.; Bogdanov, S.; Arzumanian, Z.; Gendreau, K.C.; Guillot, S.; Ho, W.C.G.; Lattimer, J.M. ; Loewenstein, M.; et al. The radius of PSR J0740+6620 from NICER and XMM-Newton data. Astrophys. J. Lett. 2021, 918, L28, 31pp.
  * (11) Salmi, T.; Vinciguerra, S.; Choudury, D.; Riley, T.E.; Watts, A.L.; Remillard, R.A.; Ray, P.S.; Bogdanov, S.; Guillot, S.; Arzumanian, Z.; et al. The radius of PSR J0740+6620 from NICER with NICER background estimates. Astrophys. J. 2022, 941, 150–173.
  * (12) Abbott, B. P.; et al. (LIGO Scientific Collaboration and Virgo Collaboration). Properties of the binary neutron star merger GW170817. Phys. Rev. 2019, X9, 011001, 32pp.
  * (13) Abbott, B. P.; et al. (LIGO Scientific Collaboration and Virgo Collaboration). GW170817: measurement of neutron star radii and equation of state. Phys. Rev. Lett. 2018, 121, 161101, 16pp.
  * (14) Abbott, B.P.; et al. (LIGO Scientific Collaboration and Virgo Collaboration). GW190425: Observation of a compact binary coalescence with total mass $\sim$ 3.4 $M_{\odot}$. Astrophys. J. Lett. 2020, 892, L3, 24pp.
  * (15) Raaijmakers, G.; Greif, S.K.; Hebeler, K.; Hinderer, T.; Nissanke, S.; Schwenk, A.; Riley, T.E.; Watts, A.; Lattimer, J.M.; Ho W.C.G. Constraints on the dense matter equation of state and neutron star properties from NICER’s mass-radius estimate of PSR J0740+6620 and multimessenger observations. Astrophys. J. Lett. 2021, 918, L29, 13pp.
  * (16) Pang, P.T.H.; Tews, I.; Coughlin, M.W.; Bulla, M.; Van Den Broeck, C.; Dietrich, T. Nuclear physics multimessenger constraints on the neutron star equation of state: adding NICER’s PSR J0740+6620 measurement. Astrophys. J. 2021, 922, 14, 9pp.
  * (17) Legred, I.; Chatziiannou, K.; Essick,R. ; Han, S.; Landry, P. Impact of the PSR J0740+6620 radius constraint on the properties of high-density matter. Phys. Rev. 2021, D104, 063003, 21pp.
  * (18) Biswas, B.; Datta, S. Constraining neutron star properties with a new equation of state insensitive approach. Phys. Rev. 2022, D106, 043012, 10pp.
  * (19) Ecker, C.; Rezzolla, L. A general, scale-independent description of the sound speed in neutron stars. Astrophys. J. Lett. 2022, 939, L35, 8pp.
  * (20) Altiparmak, S.; Ecker, C.; Rezzolla, L. On the sound speed in neutron stars. Astrophys. J. Lett. 2022, 939, L34, 9pp.
  * (21) Huth, S.; Pang, P.T.H.; Tews, I. ; Dietrich, T. ; Le Fevre, A.; Schwenk, A.; Trautmann, W.; Agarwal, K.; Bulla, M. ; Coughlin, M.W.; Van Den Broeck, C. Constraining neutron star matter with microscopic and macroscopic collisions. Nature 2022, 606, 276–280.
  * (22) Annala, E.;, Gorda, T.; Hirvonen, J.; Komoltsev, O.; Kurkela, A.; Nättilä, J.; Vuorinen, A. Strongly interacting matter exhibits deconfined behavior in massive neutron stars. 2023, arXiv:2303.11356, 12pp.
  * (23) Somasundran, R.; Tews, I.; Margueron, J. Investigating signatures of phase transitions in neutron-star cores. Phys. Rev. 2023, C107, 025801, 7pp.
  * (24) Essick, R.; Legred, I.; Chatziioannou, K.; Han, S.; Landry, P. Phase transition phenomenology with nonparametric representations of the neutron star equation of state. Phys. Rev. 2023, D108, 043013, 30pp.
  * (25) Brandes, L. ; Weise, W.; Kaiser, N. Inference of the sound speed and related properties of neutron stars. Phys. Rev. 2023, D107, 014011, 25pp.
  * (26) Lim, Y.; Holt, J.W. Neutron star radii, deformabilities, and moments of inertia from experimental and ab-initio theory constraints of the ${}^{208}Pb$ neutron skin thickness. Galaxies 2022, 10, 99, 20pp.
  * (27) Brandes, L. ; Weise, W.; Kaiser, N. Evidence against a strong first-order phase transition in neutron star cores: impact of new data. Phys. Rev. 2023, D108, 094014, 21pp.
  * (28) Mroczek, D.; Miller, M.C.; Noronha-Hostler, J.; Yunes, N. Nontrivial features in the speed of sound inside neutron stars. 2023, arXiv:2309.02345, 31pp.
  * (29) Han, M.-Z.; Huang, Y.-J.; Tang, S.-P.; Fan, Y.-Z. Plausible presence of new state in neutron stars with masses above 0.98 $M_{TOV}$. Sci. Bull. 2023, 68, 913–919.
  * (30) Drischler, C.; Han, S.; Reddy, S. Large and massive neutronn stars: implications for the sound speed within QCD of dense matter. Phys. Rev. 2022, C105, 035808, 9pp.
  * (31) Gorda, T.; Kurkela, A.; Paatelainen, R.; Säppi, S.; Vuorinen, A. Soft interactions in cold quark matter. Phys. Rev. Lett. 2021, 127, 162003,6pp.
  * (32) Komoltsev, O.; Kurkela, A. How perturbative QCD constrains the equation of state at neutron-star densities. Phys. Rev. Lett. 2022, 128202701,7pp.
  * (33) Romani, R.W.; Kandel, D. ; Filippenko, A.W.; Brink, T.G.; Zheng, W.K. PSR J0952-0607: the fastest and heaviest known galactic neutron star. Astrophys. J. Lett. 2022, 934, L17, 6pp.
  * (34) Essick, R.; Tews, I.; Landry, P.; Reddy, S.; Holz, D.E. Direct astrophysical tests of chiral effective field theory at supranuclear densities. Phys. Rev. 2020, C102, 055803, 18pp.
  * (35) Fasano, M.; Abdelsalhin, T.; Maselli, A; Ferrari, V. Constraining the neutron star equation of state using multiband independent measurements of radii and tidal deformabilities. Phys. Rev. Lett. 2019, 123, 141101,6pp.
  * (36) Akmal, A.; Pandharipande, V.R.; Ravenhall, D.G. The equation of state of nucleon matter and neutron star structure. Phys. Rev. 1998, C58, 1804-1828.
  * (37) Nättilä, J.; Miller, M.C.; Steiner, A.W.; Kajava, J.J.E.; Suleimanov, V.F.; Poutanen, J. Neutron star mass and radius measurements from atmospheric model fits to x-ray burst cooling tail spectra. Astron. & Astrophys. 2017, 608, A31, 15pp.
  * (38) Fujimoto, Y.; Fukushima, K.; McLerran, L.D.; Praszalowicz, M. Trace anomaly as signature of conformality in neutron stars. Phys. Rev. Lett. 2022, 129, 252702, 6pp.
  * (39) Marczenko, M.; McLerran, L.D.; Redlich, K.; Sasaki, C. Reaching percolation and conformal limits in neutron stars. Phys. Rev. 2023, C107, 025802, 5pp.
  * (40) Rho, M. Dense baryonic matter predicted in "pseudo-conformal model". Symmetry 2023,15, 1271, 14pp.
  * (41) Ma, Y.-L.; Yang, W.-C. Topology and emergent symmetries in dense compact star matter. Symmetry 2023,15, 776, 26pp.
  * (42) Sumiyoshi, K; Kojo, T.; Furusawa, S. Equation of state in neutron stars and supernovae. In Handbook of Nuclear Physics, Vol. 5; Tanihata, I. , Toki, H., Kajino,T., Eds.; Springer Nature: Singapore, 2023; pp. 3127–3177.
  * (43) Fukushima, K.; Hatsuda, T. The phase diagram of dense QCD. Rep. Prog. Phys. 2011, 74, 014001.
  * (44) Holt, J.W.; Rho, M.; Weise, W. Chiral symmetry and effective field theories for hadronic, nuclear and stellar matter. Phys. Reports 2016, 621, 2-75.
  * (45) McLerran, L.; Reddy, S. Quarkyonic matter and neutron stars. Phys. Rev. Lett. 2019 122, 122701, 5pp.
  * (46) Baym, G.; Furusawa, S.; Hatsuda, T.; Kojo, T.; Togashi, H. New neutron star equation of state with quark-hadron crossover. Astrophys. J.  2019, 885, 42-49.
  * (47) Fukushima, K.; Kojo, T.; Weise, W. Hard-core deconfinement and soft-surface delocalization from nuclear to quark matter. Phys. Rev. 2020, D102, 096017, 14pp.
  * (48) Kojo, T. ; Baym, G.; Hatsuda, T. Implications of NICER for neutron star matter. Astrophys. J. 2022, 934, 46-58.
  * (49) Fujimoto, Y.; Kojo, T.; McLerran, L.D. Momentum shell in quarkyonic matter from explicit duality: a solvable model analysis. 2023, arXiv:2306.04304, 6pp.
  * (50) Han, S.; Mamun, M.A.A.; Lalit, S. ; Constantinou, C.; Prakash, M. Treating quarks within neutron stars. Phys. Rev. 2019, D100, 103022, 24pp.
  * (51) Wellenhofer, C.; Holt, J.W.; Kaiser, N.; Weise, W. Nuclear thermodynamics from chiral low-momentum interactions. Phys. Rev. 2014, C89, 064009, 22pp.
  * (52) Thomas, A.W.; Weise, W. The Structure of the Nucleon; Wiley-VCH: Berlin, Germany, 2001; pp. 232–246.
  * (53) Lin, Y.H.; Hammer, H.W.; Meißner, U.-G. New insights into the nucleon’s electromagnetic structure. Phys. Rev. Lett. 2022 128, 052002, 6pp.
  * (54) Kaiser, N.; Passemar, E. Spectral functions of nucleon form factors: three-pion continua at low energies. Eur. Phys. J. 2019, A55, 16, 12pp.
  * (55) Brown, G.E.; Rho, M.; Weise, W. Phenomenological delineation of the quark-gluon structure from nucleon electromagnetic form factors. Nucl. Phys. 1986, A454, 669–690.
  * (56) Meissner, U.-G.; Kaiser, N.; Weise, W. Nucleons as Skyrme solitons with vector mesons. Nucl. Phys. 1987, A466, 685–723.
  * (57) Hill, R.J.; Kammel, P.; Marciano, W.J.; Sirlin, A. Nucleon axial radius and muonic hydrogen - a new analysis and review. Rep. Prog. Phys. 2018, 81, 096301, 36pp.
  * (58) Kharzeev, D. E. Mass radius of the proton. Phys. Rev. 2021, D104, 054015, 10pp.
  * (59) Benhar, O. Testing the paradigm of nuclear many-body theory. Particles 2023, 6, 611–621.
  * (60) Stephanov, M.; Rajagopal, K.; Shuryak, E. Signatures of the tricritical point in QCD. Phys. Rev. Lett. 1998 81, 4816–4819.
  * (61) Fukushima, K. Spectral functions in the $\sigma$ channel near the critical end point. Phys. Rev. 2003, C67, 025003, 10pp.
  * (62) Fischer, C.S.; Luecker, J.; Welzbacher, C.A. Phase structure of three and four flavor QCD. Phys. Rev. 2014, D90, 034022, 11pp.
  * (63) Asakawa, M.; Yazaki, K. Chiral restoration at finite density and temperature. Nucl. Phys. 1989, A504, 668–684.
  * (64) Klimt, S.; Lutz, M.F.M.; Weise, W. Chiral phase transition in the SU(3) Nambu and Jona-Lasinio model. Phys. Lett.1990, B249, 386–390.
  * (65) Scavenius, O.; Mocsy, A.; Mishustin, I.N.; Rischke, D.H. Chiral phase transition within effective models with constituent quarks. Phys. Rev. 2001, C64, 045202, 10pp.
  * (66) Rößner, S.; Ratti, C.; Weise, W. Polyakov loop, diquarks and the two-flavour phase diagram. Phys. Rev. 2007, D75, 034007, 15pp.
  * (67) Rößner, S.; Hell, T.; Ratti, C.; Weise, W. The chiral and deconfinement crossover transitions: PNJL model beyond mean field. Nucl. Phys. 2008, A814, 118–143.
  * (68) Hell, T.; Rößner, S.; Christoforetti, M.; Weise, W. Thermodynamics of the three-flavor nonlocal Polyakov – Nambu – Jona-Lasinio model. Phys. Rev. 2010, D81, 074023, 21pp.
  * (69) Skokov, V.; Friman, B.; Nakano, E.; Redlich, K.; Schaefer, B.-J. Vacuum fluctuations and thermodynamics od chiral models. Phys. Rev. 2010, D82, 034029, 9pp.
  * (70) Drews, M.; Weise, W. Functional renormalization group studies of nuclear and neutron matter. Prog. Part. Nucl. Phys. 2017, 93, 69–107.
  * (71) Brandes, L.; Kaiser, N.; Weise, W. Fluctuations and phases in baryonic matter. Eur. Phys. J. 2021, A57, 243, 15pp.
  * (72) Zacci, A.; Schaffner-Bielich, J. Effects of renormalizing the chiral SU(2) quark-meson model. Phys. Rev. 2018, D97, 074011, 12pp.
  * (73) Gupta, U.S.; Tiwari, V.K. Revisiting the phase structure of the Polyakov-quark-meson model in the presence of vacuum fermion fluctuation. Phys. Rev. 2012, D85, 014010, 17pp.
  * (74) Eser, J.; Blaizot, J.-P. Thermodynamics of the parity-doublet model: symmetric nuclear matter and the chiral transition. arXiv:2309.06566 2023, 37pp.
  * (75) Lonardoni, D.; Lovato, A.; Gandolfi, S.; Pederiva, F. Hyperon puzzle: hints from Quantum Monte Carlo calculations. Phys. Rev. Lett. 2015 114, 092301, 5pp.
  * (76) Gerstung, D.; Kaiser, N.; Weise, W. Hyperon-nucleon three-body forces and strangeness in neutron stars. Eur. Phys. J. 2020, A56, 175, 13pp.
  * (77) Leong, J.; Motta, T.F.; Thomas, A.W.; Guichon, P.A.M. Dense nuclear matter with phenomenological short distance repulsion. Phys. Rev. 2023, C108, 015804, 15pp.
  * (78) Baym, G.; Chin, S.A. Landau theory of relativistic Fermi liquids. Nucl. Phys. 1976, A262, 527–538.
  * (79) Friman, B.; Weise, W. Neutron star matter as a relativistic Fermi liquid. Phys. Rev. 2023, C100, 065807, 15pp.
  * (80) Baym, G.; Pethick, C. Landau Fermi-Liquid Theory: Concepts and Applications; Wiley: New York, 1991; p.117.