11institutetext: School of Physics, Southeast University, Nanjing 211189,
China 22institutetext: Department of Physics and Astronomy, Texas A$\&$M
University-Commerce, Commerce, TX 75429, USA

# Impact of Symmetry Energy on Sound Speed and Spinodal Decomposition in Dense
Neutron-Rich Matter

Nai-Bo Zhang and Bao-An Li 1122<EMAIL_ADDRESS><EMAIL_ADDRESS>
(Received: date / Revised version: date)

###### Abstract

Using a meta model for nuclear Equation of State (EOS) with its parameters
constrained by astrophysical observations and terrestrial nuclear experiments,
we examine effects of nuclear EOS especially its symmetry energy $E_{\rm
sym}(\rho)$ term on the speed of sound squared $C^{2}_{s}(\rho)$ and the
critical density $\rho_{t}$ where $C^{2}_{s}(\rho_{t})$ vanishes (indicating
the onset of spinodal decomposition) in both dense neutron-rich nucleonic
matter relevant for relativistic heavy-ion collisions and the cold $n+p+e+\mu$
matter in neutron stars at $\beta$-equilibrium. Unlike in nucleonic matter
with fixed values of the isospin asymmetry $\delta$, in neutron stars with a
density dependent isospin profile $\delta(\rho)$ determined consistently by
the $\beta$ equilibrium and charge neutrality conditions, the
$C^{2}_{s}(\rho)$ almost always show a peak and then vanishes at $\rho_{t}$.
The latter strongly depends on the high-density behavior of $E_{\rm
sym}(\rho)$ if the skewness parameter $J_{0}$ characterizing the stiffness of
high-density symmetric nuclear matter (SNM) EOS is not too far above its
currently known most probable value of about $-190$ MeV inferred from recent
Bayesian analyses of neutron star observables. Moreover, in the case of having
a super-soft $E_{\rm sym}(\rho)$ that is decreasing with increasing density
above about twice the saturation density of nuclear matter, the $\rho_{t}$ is
significantly lower than the density where the $E_{\rm sym}(\rho)$ vanishes
(indicating the onset of isospin-separation instability and pure neutron
matter formation) in neutron star cores.

###### pacs:

21.65.Mn; 26.60.Kp

##  1 Introduction

The energy per nucleon $E(\rho,\delta)$ in cold neutron-rich nucleonic matter
of density $\rho$ and isospin asymmetry $\delta=(\rho_{n}-\rho_{p})/\rho$
where $\rho_{n}$ and $\rho_{p}$ are the densities of neutrons and protons,
respectively, can be written as Bom91

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

where $E_{0}(\rho)$ is the energy per nucleon in symmetric nuclear matter
(SNM) while the symmetry energy $E_{\rm{sym}}(\rho)$ encodes the energy cost
to make the matter more neutron-rich Tesym . Both the $E_{0}(\rho)$ and
$E_{\rm{sym}}(\rho)$ are important for the stiffness and stability of neutron-
rich nuclear matter. For example, within the minimal model of neutron stars
consisting of neutrons (n), protons (p), and electrons (e) the
incompressibility of n+p+e matter at $\beta$ equilibrium is approximately
Lattimer00 ; Kubis ; JXu2
$K_{\mu}(\rho)=\rho^{2}\frac{d^{2}E_{0}}{d\rho^{2}}+2\rho\frac{dE_{0}}{d\rho}+\delta^{2}[\rho^{2}\frac{d^{2}E_{\rm
sym}}{d\rho^{2}}+2\rho\frac{dE_{\rm sym}}{d\rho}-2E^{-1}_{\rm
sym}(\rho)(\rho\frac{dE_{\rm sym}}{d\rho})^{2}].$ When $K_{\mu}(\rho)$ is
negative the corresponding speed of sound $C_{s}(\rho)$ become imaginary,
indicating the onset of the so-called spinodal decomposition. When this
happens, the density fluctuations in the dynamically unstable region grow
exponentially.

At zero temperature, the $K_{\mu}$ can become negative at both sub-saturation
and supra-saturation densities depending on the characteristics of both
$E_{0}(\rho)$ and $E_{\rm sym}(\rho)$. At low densities, the onset of spinodal
decomposition in either the n+p+e matter or purely nucleonic matter is often
used to identify the crust-core transition in neutron stars Lattimer00 ; Kubis
; JXu2  or a liquid-gas phase transition Bertsch-Siemens ; Siemens ; Lopez ;
Goodman ; Li-Siemens ; Muller  related to the nuclear multifragmentation
phenomenon in intermediate energy heavy-ion reactions Moretto ; Bauer ; Li-
Gross ; Xu93 ; Maria1 ; Maria2 ; Pawel20 . In hot dense hadronic matter formed
in relativistic heavy-ion reactions, a high-density spinodal decomposition has
been associated with a first-order hadron-quark phase transition, see, e.g.,
Refs. Risch ; Dan98 ; Dirk95 ; LiKo ; Randrup . Many potential signatures of
such transition have been proposed. For example, the resulting time delay of
expansion and/or the rapid growth of fluctuations after the spinodal
decomposition may lead to the disappearance of elliptic flow Dan98  or a
minimum in the excitation function of collective flow Dirk95 ; LiKo  or an
enhanced production of composite particles Randrup  in relativistic heavy-ion
collisions. In cold neutron star matter, when the high-density SNM EOS
$E_{0}(\rho)$ is not very stiff and the symmetry energy $E_{\rm sym}(\rho)$ is
very soft, the spinodal decomposition was also found to occur in n+p+e matter
at supra-saturation densities in the cores of neutron stars Kubis ; Li20 . It
is thus possible that there are two spinodal instability regions. While the
onset of the low-density instability indicates the crust-core transition, the
nature of the high-density one is still unclear. Kubis Kubis  speculated that
the high-density spinodal instability may indicate that the homogeneous phase
in the inner core splits into mixed two phases. The new phase may include
funny phases or solidification of the central part of stellar core. It was
further pointed out that such solid inner core may have important consequences
on the rotation of a star and its magnetic properties Kubis .

Depending on the EOS of neutron-rich matter Kubis , starting from the crust-
core transition density the $K_{\mu}$ (and the corresponding $C^{2}_{s}$ may
first increase and then decrease smoothly to zero at a critical density
$\rho_{t}$ above which the second spinodal instability happens in neutron
stars. In view of the extensive studies in the literature about the speed of
sound, here we emphasize two important distinctions of the behavior of
$C^{2}_{s}$ described above from that studied extensively earlier. Firstly,
often a sudden drop of $C^{2}_{s}$ to zero in hadronic matter is taken as an
indication of a first-order phase transition. The transition density is model
dependent and is often taken as a free parameter CSS , see, e.g., Ref. SHan
for a brief recent review. Thus, a continuously decreasing $C^{2}_{s}$ to zero
does not signal a first-order phase transition as by definition the latter is
characterized by a discontinuous, nonanalytic drop in $C^{2}_{s}$. Secondly,
there are many interesting studies on the sound speed in different phases of
dense matter or neutron stars , see, e.g., Refs. Ann19 ; Jerome22 ; Rho ;
Malik ; Loz ; Reddy ; Tan ; Dris ; Agu ; Rez ; Bra ; Kan ; Kojo ; Mot ; Larry
; Kumar ; Han22 ; Eck22 ; Yuki22  and references therein. Model predictions
and/or parameterizations of $C^{2}_{s}$ Tews18 ; Tan2  often purposely make
the $C^{2}_{s}$ reach its QCD conformal limit of 1/3 at densities far
exceeding that is reachable in neutron stars. As we shall show, depending on
the high-density EOS, the possible second spinodal instability in dense
neutron-rich matter happens long before reaching the QCD conformal limit.
Moreover, some of these studies conjectured that a peak of $C^{2}_{s}$ appears
at certain density due to some fundamentally new physics happening in the
cores of neutron stars, such as the formation of quarkyonic matter Reddy
and/or derivative contribution from the trace anomaly Yuki22 . While we do not
question these fundamentally new physics as possible origins of the peak in
$C^{2}_{s}$, within a traditionally minimal neutron star model, we show that
the $C^{2}_{s}(\rho)$ in neutron stars almost always show a peak and then
vanishes at certain $\rho_{t}$ strongly depending on the high-density behavior
of $E_{\rm sym}(\rho)$ if the skewness parameter $J_{0}$ characterizing the
stiffness of high-density SNM EOS is not too far above its currently known
most probable value of about $-190$ MeV Xie-APJ1 ; Xie-APJ2 .

While effects of nuclear EOS on the crust-core (or liquid-gas) transition
density and pressure have been extensively studied, see, e.g., reviews in Ref.
LCK08 ; Pro14 ; EPJA-review ; Pro21 , little is known about the EOS effects on
the onset density of the second spinodal instability in neither heavy-ion
reactions nor neutron stars. In this work, we study the $C^{2}_{s}(\rho)$ and
its possible second vanishing point $\rho_{t}$ in both dense neutron-rich
nucleonic matter relevant for relativistic heavy-ion collisions and the
$n+p+e+\mu$ matter in neutron stars at $\beta$-equilibrium and zero
temperature. To explore easily the whole EOS parameter space allowed by our
current knowledge, we use a meta model for nuclear EOS with its parameters
constrained by astrophysical observations and terrestrial nuclear experiments.
We focus on effects of nuclear EOS especially its symmetry energy $E_{\rm
sym}(\rho)$ term on the $C^{2}_{s}(\rho)$ and $\rho_{t}$. Because of the long
time needed for $\beta$-equilibrium due to weak interactions, it is the
partial derivative of pressure $(\partial
P(\rho,\delta)/\partial\rho)_{\delta}$ at a constant $\delta$ determines the
$C^{2}_{s}$ in relativistic heavy-ion reactions which happen on very short
time scales. While in neutron stars at $\beta-$equilibrium, it is the total
derivative $dP(\rho,\delta(\rho))/d\rho$ determines the sound speed as the
density profile $\delta(\rho)$ of the isospin asymmetry $\delta$ is already
determined by the $\beta$-equilibrium and charge neutrality conditions.

We notice that while we can find the second spinodal decomposition onset
density $\rho_{t}$ from the hadronic EOS alone and there are interesting
speculations in the literature about what may happen afterwards, in this work
we do not study neither the nature of the possible phase transition nor
properties of the new phases at higher baryon densities. Nevertheless, the
vanishing point of the speed of sound in dense hadronic matter and its
dependence on the hadronic EOS parameters are important information useful for
both nuclear physics and astrophysics.

It is also necessary to note that instabilities of nuclear matter with Skyrme
forces above saturation density have been discussed previously, see, e.g.,
Refs. Kut93 ; Mar02 ; Stone03 . Some studies consider such instabilities as
being an undesired feature of Skyrme forces in calculating the EOS within the
Skyrme-Hartree-Fock approach, see, e.g., Ref. Mar02 ; Stone03 , mostly with
respect to the desire of supporting neutron stars of masses above about 2M⊙
with only the hadronic EOS (plus the EOS of leptons to make the neutron star
matter charge neutral and stay at $\beta-$equilibrium) without considering a
possible phase transition and including the EOS of the new phase. To our best
knowledge, while it may be undesirable in the sense mentioned above, no one
has concluded that the hadronic EOS having high-density instabilities is
unphysical. In our opinion, the instability just indicates the onset of a new
phase requiring additional physics to construct the high-density EOS above the
instability onset point. Only with the latter, one can use the astrophysical
requirement of supporting neutron stars at least as heavy as the latest most
massive neutron star observed. Otherwise, the latter requirement will preclude
the EOS in neutron stars from going beyond nucleonic degrees of freedom and
the appearance of possible new phases of dense matter. Similarly, it is well
known qualitatively that instabilities driven by the onset of Delta-baryons
Lav ; Radu , appearances of hyperons and/or other particles may soften the EOS
to levels that can not support 2M⊙ neutron stars. Quantitatively, however, the
softening of the EOS depends strongly on the poorly known interactions of
these new particles with mesons, nucleons and among themselves Jiang ; Cai .
Nevertheless, the inability of supporting 2M⊙ neutron stars with the softened
EOSs does not necessarily mean the EOS used to find out when/where the
instability may occur is wrong, or new phases/particles should not appear. It
just calls for new physics to be included once these new phases/particles are
produced.

In the same spirit, in the present study we do not impose the condition that
the hadronic EOS itself has to be stiff enough to support 2M⊙ neutron stars,
although most of the EOS parameters we use are actually within the 68%
confidence boundaries of the posterior probability distribution functions
inferred from Bayesian analyses of neutron star observables under the
condition that the maximum neutron star mass to be supported by each EOS is a
Gaussian function centered at $(2.01\pm 0.04)$ M⊙ Xie-APJ1 ; Xie-APJ2 .
Moreover, we use the minimum model for neutron stars consisting of neutrons,
protons, electrons and muons without considering baryon resonances and
hyperons. As mentioned earlier, we investigate the behavior of the speed of
sound and when/if it reaches zero. Once the spinodal decomposition happens,
some new physics has to be incorporated to support 2M⊙ neutron stars. For
example, consistent with several other studies in the literature, see, e.g.,
Refs. Ang ; Tang21 ; Ang1 ; Jerome-c , using the same meta model EOS used here
it was found that the most probable hadron-quark transition density is
$1.6^{+1.2}_{-0.4}\rho_{0}$ in neutron stars from a Bayesian analysis of
available neutron star observables under the condition that the EOS has to be
stiff enough to support 2M⊙ neutron stars Xie3 . Such a low hadron-quark
transition density seems undesirable as no signature of such low transition
density has been seen in heavy-ion collisions, and one may suspect that the
low transition density may pose difficulties to support massive neutron stars.
But the Bayesian analysis indicates that all one needs is a very stiff EOS of
quark matter (namely, high speed of sound in quark matter) and a relatively
large skewness parameter $J_{0}$ of symmetric nuclear matter. Unfortunately,
the posterior probability distributions of all nucleonic EOS parameters with
or without considering the hadron-quark phase transition in neutron stars are
all consistent with our current knowledge from both terrestrial experiments
and astrophysical observations, indicating a hadron-quark duality of the
EOS/composition in understanding the currently available neutron star
observables Xie3 . In this regard, it is interesting to note that under the
assumption that nucleons are the only constituents of neutron star cores, a
very recnt Bayesian inference indicates that the current observations are
fully compatible with the nucleonic hypothesis Thi . Thus the readers should
be advised that all our results and discussions in the following should be
understood with the cautions mentioned above.

The rest of this paper is organized as follows. In the next section, we
summarize the major ingredients and available constraints of a meta EOS model
for nucleonic matter and neutron star matter. We also recall definitions of a
few relevant physical quantities. We present and discuss our results in
section 3. We then conclude in section 4.

##  2 Summary of a meta model EOS for neutron-rich matter

Meta-modeling of nuclear EOS has been widely used in both astrophysics and
nuclear physics, see, e.g., Refs. Margueron:2017eqc ; Jerome18 ;
Mondal:2022cva  for detailed discussions on the advantages and disadvantages
of such approach. In our meta model for neutron-rich matter NBZ18 ;
Zhang:2018bwq , the SNM EOS $E_{0}(\rho)$ and symmetry energy
$E_{\rm{sym}}(\rho)$ in Eq. (1) are parameterized respectively according to

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

where $E_{0}(\rho_{0})$=-15.9 MeV and

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

Such models have been widely used in Bayesian inferences of nuclear EOS from
neutron star observables or forward modelings of the structures of neutron
stars and nuclei as well as their mergers and collisions. The parameters
involved have the asymptotic boundary conditions that they become Taylor
expansion coefficients when the above parameterizations are used near the
saturation density $\rho_{0}$ of SNM and $\delta=0$. Thus, the fiducial values
of the incompressibility
$K_{0}=9\rho_{0}^{2}[\partial^{2}E_{0}(\rho)/\partial\rho^{2}]|_{\rho=\rho_{0}}$
and skewness
$J_{0}=27\rho_{0}^{3}[\partial^{3}E_{0}(\rho)/\partial\rho^{3}]|_{\rho=\rho_{0}}$
of SNM EOS $E_{0}(\rho)$ at $\rho_{0}$, the magnitude
$E_{\rm{sym}}(\rho_{0})$, slope $L=3\rho_{0}[\partial
E_{\rm{sym}}(\rho)/\partial\rho]|_{\rho=\rho_{0}}$, curvature
$K_{\rm{sym}}=9\rho_{0}^{2}[\partial^{2}E_{\rm{sym}}(\rho)/\partial\rho^{2}]|_{\rho=\rho_{0}}$
and skewness
$J_{\rm{sym}}=27\rho_{0}^{3}[\partial^{3}E_{\rm{sym}}(\rho)/\partial\rho^{3}]|_{\rho=\rho_{0}}$
of nuclear symmetry energy $E_{\rm sym}(\rho)$ at $\rho_{0}$, respectively,
from astrophysical observations and terrestrial nuclear experiments provide a
guide and some constraints on the parameters of the meta model EOS.

Figure 1: An updated compilation of the slope paramter L of nuclear symmetry
energy extracted from analyses of several recent terrestrial experimenst in
comparison with the 2016 fiducial value of $L=58.7\pm 28$ MeV (light mint
band) as5  and the average value of $L=58\pm 19$ MeV (black horizonal line)
from 24 analyses (listed in Ref. LiBA21 ) of neutron star observables after
GW170817. See the text for more details.

For studying the high-density spinodal decomposition in this work, we fix the
low-density EOS parameters at their currently known most probably values,
e.g., $K_{0}=240$ MeV Garg18  and $E_{\rm{sym}}(\rho_{0})=31.7$ MeV Li2013
and $L=60\pm 20$ MeV, but vary the high-density parameters within their
current uncertain ranges if known. To justify the $L$ value we use, in Fig. 1
we provide an updated compilation of $L$ with respect to its 2016 fiducial
value of $L=58.7\pm 28$ MeV (light mint band) from averaging results of 53
independent analyses of various terrestrial and astrophysical data available
as5 . Starting from the left are the $L$ values from (1) 24 independent
analyses of neutron star observables carried out by various groups between
2017 and 2021, they give an average $L=58\pm 19$ MeV (horizontal line) LiBA21
; (2) the original analysis of the PREX-II data prex2  by Reed et al. Reed21
and 3 independent analyses of PREX-II data together with different
combinations of terrestrial and/or astrophysical data by Reinhard et al.
Reinhard1 , Essick et al. Essick21  and Yue et al. Yue22 , respectively; (3)
liquid drop model analyses using separately the PREX-I & II data together,
CREX data only, and the combination of all data assuming they are equally
reliable by Lattimer Lattimer23 ; (4) charged pion ration in heavy-ion
reactions at RIKEN by Estee et al. Estee21 ; (5) the ratio of average
transverse momentum (sky blue star) and the ratio of charged particle
multiplicites (black star) in isobar collisions (96Zr+96Zr and 96Ru+96Ru) from
STAR/RHIC experiments analysed by Xu et al. star ; (6) the difference of
charge radii of the mirror pair 54Ni-54Fe by Pineda et al. mirror ; the chiral
EFT prediction by Dirschler et al. ceft . It is beyond our ability to comment
on the underlying physics differences of the different analyses listed above
and the compilation is certainly incomplete as we only used selected results
from analyzing some new terrestrial experiments over the last 2 years.
However, we notice that the great interest stimulated by the PREX-II result
has now been increased even further by the recent publication of the CREX
result crex . Of course, there are diverse opinions in understanding the PREX
and CREX results in the community. For example, on one hand, Ref. Reinhard2
cautioned that “Until the tension between theory and experiment, or between
the two measurements, is resolved, one should exercise extreme caution when
interpreting the new APV measurements in the context of neutron skins or
nuclear symmetry energy”. On the other hand, it was found that “The two
experiments separately predict incompatible ranges of $L$ ($L=-5\pm 40$ MeV
from CREX and $L=121\pm 47$ MeV from combined PREXI & II data, respectively),
but accepting both measurements to be equally valid suggests $L=50\pm 12$ MeV,
nearly the range suggested by either nuclear mass measurements or neutron
matter theory, and is also consistent with nuclear dipole polarizability
measurements” Lattimer23 . There are also other efforts to simultaneously
explain both PREX-II and CREX data within the same theory/model with very
limited success, see, e.g., Refs. Mondal:2022cva ; Zhang22 ; Tag22 ; Yuk22 ;
Bis21 ; Cor22 . Given the results described above and shown in Fig. 1, it is
seen that the 2016 fiducial value, the average L from the 24 analyses of
neutron stars after GW179817 as well as several new results from the last 2
years are all in perfect agreement with the chiral EFT prediction of
$L=59.8\pm 4.1$ MeV ceft , albeit with different error bars. Adopting the
philosophy that all published experimental data are eqully reliable within
their reported confidence boundaries, the value of $L=60\pm 20$ MeV we use in
this work is completely consistent with our current best knowledge. Similarly,
we notice that the magnitude of symmetry energy $E_{\rm sym}(\rho_{0})=31.6\pm
2.7$ MeV Li2013  or $E_{\rm sym}(\rho_{0})=31.7\pm 3.2$ MeV as5  based on
earlier analyses of mostly terrestrial experiments also agrees very well with
the chiral EFT result ($E_{\rm sym}(\rho_{0})=31.7\pm 1.1$ MeV) ceft . These
fiducial values have been used already in many of the analyses in extracting
the $L$ values.

While several analyses have found strong supports for $K_{\mathrm{sym}}$ to be
around $-100\pm 100$ MeV Jerome22 ; Jerome18 ; LiBA21 ; Mondal , a larger
range is still often used in the literature. As to the $J_{\rm{sym}}$
parameter, only some theoretical predictions put it in the range between about
-200 to +800 MeV Dut12 ; Dutra2014 ; Tews17 ; Zhang17 . To our best knowledge,
none of the analyses of neither astrophysical nor nuclear physics data has
provided a reliable constraint on $J_{\rm{sym}}$. Consequently, the high-
density behavior of $E_{\rm sym}(\rho)$ above about $2\rho_{0}$ is largely
open LiBA21 ; Zhang21 .

To see how diverse the metal model $E_{\rm{sym}}(\rho)$ can be, a few examples
are shown in Fig. 2 by varying the $L$, $K_{\mathrm{sym}}$, and
$J_{\mathrm{sym}}$ parameters in their uncertainty ranges discussed above.

Figure 2: The density dependences of nuclear symmetry energy with different
$L$, $K_{\rm sym}$, and $J_{\rm sym}$ values as indicated.

It is seen that the generated $E_{\rm{sym}}(\rho)$ ranges from being super-
soft (first increases up to about $2\rho_{0}$ then start decreasing with
increasing density and can even become negative at high densities) to super-
stiff (increases with increasing density much fast than linear). The spread of
nuclear symmetry energy at suprasaturation densities is consistent with the
present constraints Zhang:2018bwq  and predictions from different models EPJA-
review ; LiBA21 ; Chen18 . Interestingly, one can even simulate a “cusp” in
$E_{\rm sym}(\rho)$ due to a topology change as predicted by a pseudo-
conformal symmetry theory Rho . To our best knowledge, currently neither
experimental findings nor fundamental physics principles can rule out any of
the high-density behaviors of $E_{\rm sym}(\rho)$ illustrated in Fig. 2. It
has been recognized that the poorly known isospin-dependent tensor force and
the corresponding short-range nucleon-nucleon correlation as well as the spin-
isospin dependence of three-body force in neutron-rich matter are among the
most important origins of the wide open behavior of high-density symmetry
energy EPJA-review ; Cai22 . Impacts of the diverse high-density behaviors of
$E_{\rm sym}(\rho)$ on properties of neutron stars and heavy-ion collisions
have also been studied in several works, see, e.g., Refs. Rho ; Kut93 ; LiBA21
; Cai22 ; Kut94 ; LiBA02 ; Kut93a ; Szm06 ; wen ; Kho96 ; Bas07 ; Ban00 ;
Kubis1 . Is there any physical meaning for the spread in the density
dependence of symmetry energy beyond about twice the nuclear saturation
density? Yes, albeit there is currently no consensus and model-independent
answer, see, e.g., Ref. EPJA-review  for a recent review. In fact, to pin down
the physics origin of the spreading of symmetry energy at suprasaturation
densities and find sensitive observables to probe it has been a longstanding
goal of many theoretical and experimental efforts. The latest and most
detailed discussions on this question can be found in subsection VI-D of a
very recent white paper on “Dense Nuclear Matter Equation of State from Heavy-
Ion Collisions” Sorensen .

On the other hand, there are several available constraints on the $J_{0}$
parameter which characterizes the stiffness of SNM EOS at densities above
about $(2-3)\rho_{0}$. For example, recent Bayesian analyses of neutron star
radii from LIGO and NICER observations have inferred the most probable value
of $J_{0}$ to be about $J_{0}=-190_{-40}^{+40}$ at 68% confidence level Xie-
APJ1 ; Xie-APJ2 . While its most probable value extracted from a Bayesian
analysis of nuclear collective flow in relativistic heavy-ion collisions is
$J_{0}=-180^{+100}_{-110}$ MeV at 68% confidence level XieLi-JPG . Moreover,
it was found that to support an earlier reported neutron star of mass 2.17 M⊙
a minimum value of about $J_{0}=-180$ MeV is needed NBZ217 . While the
causality requirement can provide an upper limit on $J_{0}$, the range of
about $-90<J_{0}<200$ MeV obtained in this way depends strongly on the
correlation between the high-density SNM EOS and symmetry energy NBZ217 . We
will thus explore effects of $J_{0}$ in a large range around its most probable
value of about -190 MeV.

It is worth noting here that most of the microscopic models/theories have a
more complex behavior with density compared to the meta model EOS used here.
In particular, higher orders in powers of the density and, therefore, the same
coefficients may have different meanings in different models/theories. For
example, most Skyrme models predict $J_{0}<-300$ MeV, some RMF models may
predict a very negative value as the Skryme forces, or large positive values,
and ab-initio theories predict various negative values. Concerning the $J_{\rm
sym}$ all models predict quite large values except for chiral EFT, see, e.g.,
Refs. Margueron:2017eqc ; Duc  for more detailed discussions. We are not
knowledgeable enough to explain why the $J_{0}$ and $J_{\rm sym}$ from these
different models/theories may be very different. In fact, to understand and
constrain their predictions in a common framework has been a main science
driver in the field. We notice that (1) each one of these models/theories has
multiple parameters and they are often strongly correlated, (2) often the
structures/formalisms of these models/theories are not flexible enough to be
used directly in Bayesian-kind of statistical inferences, (3) one of the main
advantages of meta model EOSs is to fully explore the whole EOS parameter
space freely without the structural limitations that most microscopic many-
body models/theories have, (4) the most probable values and uncertainty ranges
of $J_{0}$ and $J_{\rm sym}$ discussed above are from both Bayesian inferences
and bruteforce inversions of neutron star observables using the same meta EOS
model used here. Thus, we have at least self-consistency although we can not
explain why some other EOS models/theories may give different values for
$J_{0}$ and $J_{\rm sym}$. Considering this situation, we also show as
demonstrations in several cases results from using $J_{0}$ and $J_{\rm sym}$
values significantly different from their most probable values given above.

The energy density in neutron stars consisting of neutrons, protons,
electrons, and muons at $\beta$-equilibrium is given by

$\epsilon(\rho,\delta)=\rho[E(\rho,\delta)+M_{N}]+\epsilon_{l}(\rho,\delta),$
(4)

where $M_{N}$ represents the average nucleon mass and
$\epsilon_{l}(\rho,\delta)$ denotes the lepton energy density. The pressure is
then

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

The pressure in neutron stars at $\beta$-equilibrium becomes barotropic once
the density-profile of the isospin asymmetry $\delta(\rho)$ is obtained self-
consistently from the $\beta$-equilibrium condition

$\mu_{n}-\mu_{p}=\mu_{e}=\mu_{\mu}\approx 4\delta E_{\rm{sym}}(\rho)$ (6)

and the charge neutrality condition $\rho_{p}=\rho_{e}+\rho_{\mu}.$ The
chemical potential $\mu_{i}$ for a particle $i$ is obtained from
$\mu_{i}=\partial\epsilon(\rho,\delta)/\partial\rho_{i}$ Wiringa88 ;
Lattimer91 . Note that the last term in Eq. (6) is obtained assuming the mass
of electrons can be neglected and the fraction of muons is smaller than that
of electrons. We notice that in nucleonic matter at fixed isospin asymmetries
relevant for heavy-ion reactions, the nucleonic part of the above equations
gives the corresponding pressure $P(\rho,\delta)$.

For easy of our following discussions, it is useful to recall that the
pressure in the n+p+e model of neutron stars is given by Lat01

$P(\rho,\delta)=\rho^{2}[\frac{dE_{\rm{SNM}}(\rho)}{d\rho}+\frac{dE_{\rm{sym}}(\rho)}{d\rho}\delta^{2}]+\frac{1}{2}\delta(1-\delta)\rho
E_{\rm sym}(\rho)$ (7)

and the isospin profile $\delta(\rho)$ (or the corresponding proton fraction
$x_{p}(\rho)=(1-\delta)/2$) at $\beta$-equilibrium is determined completely by
the $E_{\rm sym}(\rho)$ via Lat01

$\displaystyle x_{p}(\rho)=0.048\left[E_{\rm sym}(\rho)/E_{\rm
sym}(\rho_{0})\right]^{3}(\rho/\rho_{0})(1-2x_{p}(\rho))^{3}.$ (8)

Based on this equation and as we shall illustrate numerically, as the $E_{\rm
sym}(\rho)$ goes to zero at high densities in the case of a super-soft
symmetry energy, the $\delta$ approaches 1 and if the $E_{\rm sym}(\rho)$ is
super-stiff the resulting $\delta$ will vanish at very high densities.

Figure 3: The speed of sound $C^{2}_{NM}$ in unit of $c^{2}$ in nucleonic
matter with fixed isospin asymmetries as a function of density using the
symmetry energy $E_{\rm sym}(\rho)$ functions with $J_{\rm{sym}}$ varying
between -200 and 800 MeV as shown in the right panel of Fig. 2. The orange
dashed line corresponds to the conformal limit $C^{2}<1/3$.

As the speed of sound $C_{s}(\rho)$ is a measure of the stiffness of nuclear
EOS, its high-density behavior is largely controlled by the skewness $J_{0}$
of SNM EOS as well as the curvature $K_{\mathrm{sym}}$ and skewness
$J_{\rm{sym}}$ of nuclear symmetry energy. Through the charge neutrality and
$\beta$-equilibrium conditions, the $K_{\mathrm{sym}}$ and $J_{\rm{sym}}$ also
determine the high-density behavior of the isospin profile $\delta(\rho)$ of
neutron stars. As we shall discuss, the $\delta(\rho)$ plays a significant
role in determining the density dependence of the sound speed $C_{s}(\rho)$
and its possible vanishing point $\rho_{t}$. In the following, we denote the
speed of sound squared in nucleonic matter of fixed isospin asymmetry $\delta$
with $C^{2}_{NM}$ (the so-called adiabatic speed), and that in neutron stars
at $\beta$-equilibrium with $C^{2}_{s}$ (the so-called equilibrium sound
speed). More specifically, they are respectively defined by

$C^{2}_{NM}\equiv\left({\partial
P\over\partial\epsilon}\right)_{\delta}~{}{\rm
and}~{}~{}C^{2}_{s}\equiv{dP\over d\epsilon}.$ (9)

It has been shown earlier that in terms of the pressure $P$ and
incompressibility $K\equiv 9\left({\partial
P\over\partial\rho}\right)_{\delta}$, the adiabatic sound speed in unit of
$c^{2}$ in nuclear matter can be rewritten as Jerome22 ; Blaizot ; DongLai

$C^{2}_{NM}=\frac{K}{9(M_{N}+E(\rho,\delta)+P/\rho)}.$ (10)

Since the density profile $\delta(\rho)$ has to be calculated numerically
above the $\mu$ production threshold, no analytical expression for the
$C^{2}_{s}$ is available.

We notice that besides the equilibrium sound speed one can also calculate the
adiabatic sound speed in studying the core g-mode of neutron star oscillations
DongLai ; RG1  and its dependence on the $E_{\rm sym}(\rho)$ Malik ; Loz ;
EPJA-review ; DongLai . While a fluid element in a neutron star in a core
g-mode oscillation is always in mechanical equilibrium with the surroundings,
its composition stays a constant during the adiabatic oscillation process
because the weak interaction timescale leading to $\beta$-equilibrium is much
longer than the dynamical timescale of the oscillation. Namely, the
established density profile $\delta(\rho)$ of neutron stars at
$\beta$-equilibrium will not change with time during the g-mode oscillation.
As discussed in detail in Ref. DongLai , the adiabatic speed of sound is then
determined by the partial derivative of pressure with respect to density at a
constant proton/electron fraction or isospin asymmetry $\delta$ (the same
expression as for $C^{2}_{NM}$ but with the pressure $P$ for neutron stars).
The difference between the adiabatic and the $\beta$-equilibrium sound speeds
determines the g-mode frequency. It also determines the so-called Brunt-
Väisälä frequency quantifying the stability condition of the core g-mode
against convective instability DongLai ; RG1 . Interestingly, the high-density
behavior of $E_{\rm sym}(\rho)$ was found to have a significant effect on the
core g-mode frequency, see, e.g., Ref. DongLai  and the review in Section 7.1
of Ref. EPJA-review .

As noticed earlier, there are many interesting studies on the sound speed in
different phases of dense matter or neutron stars under various conditions,
see, e.g., Refs. Ann19 ; Jerome22 ; Rho ; Malik ; Loz ; Reddy ; Tan ; Dris ;
Agu ; Rez ; Bra ; Kan ; Kojo ; Mot ; Larry ; Kumar ; Han22 ; Eck22 ; Yuki22 .
As outlined above, we focus on studying the density dependence of the sound
speed in purely nucleonic matter or the $n+p+e+\mu$ matter and its possible
vanishing point $\rho_{t}$. The latter may be used as an input in constructing
more complicated EOSs to study some of the interesting and unresolved issues
about the phase diagram of dense neutron-rich matter.

##  3 Results and Discussions

In the following, we present our results on the speeds of sound in nucleonic
matter and neutron star matter separately. We examine effects of the EOS
parameters on the density dependence of the speeds of sound and their possible
vanishing points.

###  3.1 Speed of sound in neutron-rich nucleonic matter with fixed isospin
asymmetries

Shown in Fig. 3 is the speed of sound $C^{2}_{NM}$ in nucleonic matter with
fixed isospin asymmetries of $\delta=0.1,0.3$ and 0.6, respectively, as a
function of density using the symmetry energy $E_{\rm sym}(\rho)$ functions
with $J_{\rm{sym}}$ varying between -200 and 800 MeV as shown in the right
panel of Fig. 2. All other EOS parameters are fixed at their currently known
most probable values as we discussed earlier. As one expects, as the matter
becomes more neutron-rich, effects of $J_{\rm{sym}}$ become larger.
Interestingly, with relatively small $\delta$ and/or $J_{\rm{sym}}$ values the
$C^{2}_{NM}$ shows a peak. From inspecting the expression of Eq. (10) for
$C^{2}_{NM}$ in the limiting case of SNM, one can see that while the
denominator keeps increasing with density, the incompressibility $K$ in the
numerator peaks at certain density since the $K_{0}$ is positive but $J_{0}$
is negative. Thus, the $C^{2}_{NM}$ in SNM normally has a peak. With the
parameters $L$ and $K_{\rm{sym}}$ fixed, the isospin asymmetric energy,
pressure, and its derivative all increase with positive and increasing
$J_{\rm{sym}}$. Their net effect is adding a positive and continuously
increasing contribution to $C^{2}_{NM}$, making its peak disappears gradually
when the $E_{\rm sym}(\rho)$ becomes very stiff for large $J_{\rm{sym}}$
values. Of course, if one uses negative $J_{\rm{sym}}$ values leading to
super-soft $E_{\rm sym}(\rho)$ functions, the $C^{2}_{NM}$ not only has a peak
but also vanishes at a lower density $\rho_{t}$ as we shall discuss next.

Figure 4: A comparison of the transition density $\rho_{t}$ where
$C^{2}_{NM}(\rho_{t})$=0 and the critical density $\rho_{c}$ where $E_{\rm
sym}(\rho_{c})=0$ indicating the onset of isospin separation instability with
EOS parameters leading to super-soft symmetry energies.

As shown in Fig. 2, with small or negative $J_{\rm{sym}}$ values, the $E_{\rm
sym}(\rho)$ is super-soft and vanishes at some point or becomes negative at
higher densities. It is thus interesting to compare the transition density
$\rho_{t}$ where $C^{2}_{NM}(\rho_{t})$=0 and the critical density $\rho_{c}$
where $E_{\rm sym}(\rho_{c})=0$. It is known that a vanishing $E_{\rm
sym}(\rho)$ in nucleonic matter indicates the onset of the so-called isospin
separation instability, namely it is energetically more favorable to split SNM
into pure neutron matter and proton matter when the $E_{\rm sym}(\rho)$
becomes negative Kut94 ; LiBA02 . For neutron stars, two isospin separation
mechanisms were considered by Kutschera et al. Kut94 : A bulk separatmn of
protons and neutrons or formation of a neutron bubble around a single
localized proton (the so-called proton polaron). The latter was considered to
be more likely to occur than the bulk separtion in neutron star matter.
Moreover, it was speculated that the proton polarons in a neutron star core
could form a lattice due to the Coulomb interaction at low temperature.
Properties of such matter has been studied in a series of papers, see, e.g.,
Refs. Kut93 ; Kut93a ; Szm06 . Possible ramifications of isospin separation
instability in heavy-ion reactions and neutron stars have also been studied in
several other works, see., e.g., Refs. Rho ; LiBA21 ; Cai22 ; LiBA02 ; wen ;
Kho96 ; Bas07 ; Ban00 ; Kubis1 . However, many interesting issues remain to be
addressed.

Figure 5: Effects of $J_{0}$ characterizing the stiffness of SNM EOS on the
speed of sound in nucleonic matter with the $E_{\rm sym}(\rho)$ parameters
fixed at the values specified leading to a stiff $E_{\rm sym}(\rho)$ function.
The solid, dashed and dot-dashed lines are for $\delta=0.1,0.3$ and 0.6,
respectively. The orange dashed line corresponds to the conformal limit
$C^{2}<1/3$.

Shown in Fig. 4 is a comparison of the two transition densities where either
$C^{2}_{NM}(\rho_{t})$=0 or $E_{\rm sym}(\rho_{c})=0$ as a function of
$K_{\rm{sym}}$ with several small $J_{\rm{sym}}$ values all leading to super-
soft $E_{\rm sym}(\rho)$ functions. As illustrated in Fig. 2, it is the
combination of $K_{\rm{sym}}$ and $J_{\rm{sym}}$ that is determining the high-
density behavior of $E_{\rm sym}(\rho)$ above about $2\rho_{0}$. It is seen
that when both the $K_{\rm{sym}}$ and $J_{\rm{sym}}$ are small leading to very
soft $E_{\rm sym}(\rho)$ functions that vanish quickly, the critical density
for $E_{\rm sym}(\rho_{c})=0$ (isospin separation instability) is smaller than
that for $C^{2}_{NM}(\rho_{t})$=0 (spinodal decomposition instability). Since
most of the existing analyses support $K_{\rm{sym}}\approx-100\pm 100$ MeV at
68% confidence level and if indeed the $J_{\rm{sym}}$ is small, a super-soft
$E_{\rm sym}(\rho)$ is realized, and then the isospin separation occurs first.
However, as we discussed earlier, there is currently no experimental
constraint on the $J_{\rm{sym}}$, leaving all possibilities open. Thus, our
findings here add more importance for constraining the high-density behavior
of nuclear symmetry energy especially its skewness parameter $J_{\rm{sym}}$.

In the results presented above, we have used the currently known most probable
value of $-190$ MeV for the skewness parameter $J_{0}$ of SNM EOS. As we
discussed earlier, the $J_{0}$ still has a large uncertainty range. To reveal
its quantitative effects on the speed of sound, we vary $J_{0}$ between 0 and
-190 MeV using a relatively stiff $E_{\rm sym}(\rho)$ function created with
$K_{\rm{sym}}=-230$ MeV and $J_{\rm{sym}}=800$ MeV. The results are shown in
Fig. 5. It is seen that the peak in the speed of sound gradually disappears if
the $J_{0}$ is far from its most probable value of -190 MeV such that the SNM
EOS becomes very stiff at high densities. This effect is stronger for more
neutron-rich matter as a very large $J_{\rm{sym}}$ of 800 MeV is used. These
features are easily understood from the terms involved in the expression for
$C^{2}_{NM}$ in Eq. (10).

In short, in nucleonic matter with fixed isospin asymmetries, if the $J_{0}$
is close to its most probable value of -190 MeV and the $E_{\rm sym}(\rho)$ is
not too stiff, the $C^{2}_{NM}(\rho)$ shows a peak and may become zero at high
densities. However, if the $E_{\rm sym}(\rho)$ is super-stiff especially with
high isospin asymmetries the $C^{2}_{NM}(\rho)$ will keep increasing with
density until reaching the causality limit. As we shall illustrate and discuss
next, the physical consequence of the last case will be changed completely
when the $\beta$-equilibrium condition is introduced in building the EOS for
neutron stars.

###  3.2 Speed of sound in charge neutral $npe\mu$ matter in neutron stars at
$\beta$-equilibrium

Besides the contributions from leptons, the most important difference between
the EOSs of nucleonic matter at fixed $\delta$ and neutron star matter is the
consequence of charge neutrality and $\beta$-equilibrium in the latter. It is
well known that in neutron stars at $\beta$-equilibrium, the isospin profile
$\delta(\rho)$ is uniquely determined by the Eq.(6). As mentioned earlier, the
resulting $\delta(\rho)$ approaches one (zero) when the $E_{\rm sym}(\rho)$
becomes zero (very high or stiff). This result is also easy to understand
qualitatively from the energy conservation point of view. As shown in Eq. (2),
the $E(\rho,\delta)\propto E_{\rm{sym}}(\rho)\cdot\delta^{2}$. To conserve the
total energy, if the $E_{\rm{sym}}(\rho)$ is high (low) at certain density
$\rho$, the $\delta(\rho)$ there will be low (high). In both cases, the same
isospin asymmetry $\delta$ of the whole system is just being distributed to
different density regions according to the well-known isospin fractionation
mechanism in asymmetric nuclear matter Muller ; liko ; baran ; shi ; HXu .
Considering two connected regions with local densities $\rho_{1}$ and
$\rho_{2}$, respectively, the chemical equilibrium condition requires that shi
$E_{\rm sym}(\rho_{1})\delta(\rho_{1})=E_{\rm sym}(\rho_{2})\delta(\rho_{2})$.
Thus, for a given $E_{\rm sym}(\rho)$ function, the local isospin asymmetries
will adjust themselves according to the relative symmetry energy in those
regions.

Figure 6: The density profile of isospin asymmetry $\delta(\rho)$ (upper
panel) and the corresponding sound speed squared $C^{2}_{s}(\rho)$ in unit of
$c^{2}$ in neutron stars at $\beta$-equilibrium with $J_{\rm{sym}}$ varying
between -200 and 800 MeV but other parameters fixed at their currently known
most probable values indicated. The orange dashed line corresponds to the
conformal limit $C^{2}<1/3$.

Shown in Fig. 6 are the density profile of isospin asymmetry $\delta(\rho)$
(upper panel) and the corresponding squared speed of sound in neutron stars at
$\beta$-equilibrium with $J_{\rm{sym}}$ varying between -200 and 800 MeV but
other parameters fixed at their currently known most probable values. Again,
the corresponding $E_{\rm sym}(\rho)$ is the one shown in the right panel of
Fig. 2. Indeed, it is seen that the super-stiff $E_{\rm sym}(\rho)$ with
$J_{\rm{sym}}$=800 MeV leads to a vanishing $\delta$ (SNM) at high densities.
While the super-soft $E_{\rm sym}(\rho)$ with $J_{\rm{sym}}\leq$0 leads
quickly to $\delta=1$ (pure neutron matter) at high densities. Namely, in both
cases the product of $E_{\rm sym}(\rho)$ and $\delta$ at high densities reach
zero as required by the $\beta$-equilibrium and charge neutrality conditions.
Consequently, regardless of the high-density behavior of $E_{\rm sym}(\rho)$
(super-soft or super-stiff), the equilibrium speed of sound $C^{2}_{s}(\rho)$
in neutron stars at $\beta$-equilibrium always show a peak with its position
depends on the values of $J_{\rm{sym}}$ and other EOS parameters as we shall
discuss in more detail. Therefore, compared to the results shown in the
previous subsection, it is clearly seen that the high-density behaviors of
sound speeds in neutron stars at $\beta$-equilibrium and nucleonic matter with
fixed $\delta$ values are very different.

Figure 7: Upper panel: A comparison of the critical ($\rho_{c}/\rho_{0}$) and
thansition ($\rho_{t}/\rho_{0}$) densities where the symmetry energy $E_{\rm
sym}(\rho_{c})=0$ with that where the sound speed $C^{2}_{s}(\rho_{t})=0$ in
neutron stars at $\beta-$equilibrium as functions of $K_{\rm sym}$ with
$J_{\rm sym}=0,-100$ and -200 MeV, respectively. Lower panel: The transition
asymmetry $\delta_{t}$ at $\rho_{t}$ where $C^{2}_{s}(\rho_{t})=0$ as
functions of $K_{\rm sym}$ with $J_{\rm sym}=0,-100$ and -200 MeV. The
horizontal dotted line corresponds to the critical $\delta$ below which the
direct URCA process happens. All other parameters are fixed at their currently
known most probable values.

Again, it is interesting to compare the critical density $\rho_{c}/\rho_{0}$
where the symmetry energy $E_{\rm sym}(\rho_{c})=0$ with $\rho_{t}/\rho_{0}$
where the sound speed $C^{2}_{s}(\rho_{t})=0$ in the case of having super-soft
$E_{\rm sym}(\rho)$ functions. Shown in the upper panel of Fig. 7 are the two
critical densities as functions of $K_{\rm sym}$ with $J_{\rm sym}=0,-100$ and
-200 MeV, respectively, and all other parameters fixed at their currently
known most probable values. The horizontal dotted line corresponds to the
critical $\delta_{t}$ below which the direct URCA process happens. Several
interesting observations and speculations can be made: (1) since the
$\delta_{t}$ becomes less than the direct URCA limit mostly with
large/positive $K_{\rm sym}$ values far above its most probable value, the
direct URCA is unlikely to happen before the spinodal decomposition sets in.
To our best knowledge, currently there is no strong observational evidence for
direct URCA and most of the cooling curves can be well described by the slow
modified URCA processes within normally very large error bars. We thus can not
say much about the implication of our work here on cooling of neutron stars,
(2) in the whole EOS parameter space considered, the spinodal decomposition
onset density $\rho_{t}$ is always lower than the critical density $\rho_{c}$
where $E_{\rm sym}(\rho_{c})=0$ when the isospin separation instability sets
in. Around the most probable value of $K_{\rm sym}=-100$ MeV, as seen in the
upper panel of Fig. 7, the spinodal decomposition density is around the outer
core density while the isospin separation density is around the inner core
density. Thus, we speculate that these two instabilities may happen in the
outer and inner core regions of neutron stars separately. Consequences and/or
observational evidences of this speculation deserves further studies. We
emphasize again that the isospin separation instability only happens when the
high-density symmetry energy is super-soft. On the other hand, regardless of
the high-density behavior (super-soft or super-stiff) of symmetry energy, the
speed of sound almost always first goes up, then goes down and finally
vanishes at $\rho_{t}$ in neutron stars. Of course, as we emphasized earlier,
at this time we do not know what new phases may appear and how to handle them
above $\rho_{t}$.

Figure 8: The spinodal decomposition onset density $\rho_{t}$ in 3D in neutron
stars at $\beta-$equilibrium as a function of $J_{\rm sym}$ and $K_{\rm sym}$
with $L=$ 40, 60, 80 and 140 MeV, respectively. All other parameters are fixed
at their currently known most probable values.

Figure 9: Effects of the skewness $J_{0}$ of SNM EOS on the spinodal
decomposition onset density $\rho_{t}$ in neutron stars at
$\beta$-equilibrium.

###  3.3 Onset density $\rho_{t}$ of spinodal decomposition in neutron stars
at $\beta$-equilibrium

We now focus on examining the spinodal decomposition onset density $\rho_{t}$
where $C^{2}_{s}(\rho_{t})=0$ in neutron stars at $\beta$-equilibrium. First,
we fix the $J_{0}$ at its most probable value of -190 MeV and study the
dependence of $\rho_{t}$ on the slope $L$, curvature $K_{\rm sym}$, and
skewness $J_{\rm sym}$ parameters of $E_{\rm sym}(\rho)$ within their current
uncertainty ranges in Fig. 8. Overall, a softer $E_{\rm sym}(\rho)$ leads to a
smaller $\rho_{t}$. Approaching the lower limits of $K_{\rm sym}$ and $J_{\rm
sym}$ where the $E_{\rm sym}(\rho)$ is super-soft, the $\rho_{t}$ can be as
low as about $2\rho_{0}$ depending on the L value used. Interestingly, such a
low onset density of spinodal decomposition is actually consistent with the
lower boundary of hadron-quark transition density inferred from several recent
Bayesian analyses of neutron star observables Ang ; Tang21 ; Jerome-c ; Xie3 ;
Ferr . However, at the other limit where both the $K_{\rm sym}$ and $J_{\rm
sym}$ take large positive values, the $E_{\rm sym}(\rho)$ is super-stiff and
the resulting $\rho_{t}$ are larger than $6\rho_{0}$. It is also seen that
varying the $L$ parameter from 40 to 80 MeV causes an approximately 30% change
in $\rho_{t}$. The PREX-II result indicates $L=106\pm 37$ MeV Reed21 . While
its mean is far above the fiducial value of L as shown in Fig. 1, its lower
boundary is in the 2016 fiducial band of L. To see implications of its upper
boundary on the transition density $\rho_{t}$, results with $L=140$ MeV is
also shown in Fig. 8. Obviously, it makes the symmetry energy very stiff with
the same $K_{\rm sym}$ and $J_{\rm sym}$ parameters compared to the results
using lower L values. Consequently, the $\rho_{t}$ becomes much higher as one
expects. The wide range of $\rho_{t}$ due to the uncertainty of the $E_{\rm
sym}(\rho)$ especially at high densities signifies again the need to better
constrain the density dependence of nuclear symmetry energy.

Next, we examine in Fig. 9 effects of $J_{0}$ on $\rho_{t}$ while the $E_{\rm
sym}(\rho)$ parameters are also varied around their most probable values. For
exploration purposes only, in one of the calculations we purposely used a very
low $J_{0}$ of -290 MeV leading to an EOS that is not stiff enough to support
a neutron star of mass 2.01 M⊙. We also extended the lower limit of $K_{\rm
sym}$ to -400 MeV that is far below its $1\sigma$ lower boundary of $K_{\rm
sym}\approx-100\pm 100$ MeV from several surveys of its constraints in the
literature LiBA21 . It is seen that even with such a low $K_{\rm sym}$ (which
controls the behavior of $E_{\rm sym}(\rho)$ around $(1-2)\rho_{0}$) value,
the spinodal decomposition may still happen at densities below about
$(1.5-3)\rho_{0}$. However, as the $K_{\rm sym}$ increases in the case of
using a very stiff SNM EOS (e.g., with $J_{0}=-90$ MeV shown in the left panel
far above its most likely value of $-190$ MeV used in the middle panel) and a
large $J_{\rm sym}$, it may not occur as indicated by the ending of the
$\rho_{t}$ curves with $J_{\rm sym}$=200 MeV in the left panel.

At the other extreme of using a very soft SNM EOS with $J_{0}=-290$ MeV shown
in the right panel), the $\rho_{t}$ is again dominated by the SNM EOS and
becomes less dependent on the $E_{\rm sym}(\rho)$ as indicated by the tendency
of saturation at the high $K_{\rm sym}$ limit in the right panel in Fig. 9.
These features indicate that there are strong interplays and correlations
among the SNM EOS and $E_{\rm sym}(\rho)$ parameters as one expects.

In short, unlike in nucleonic matter with fixed $\delta$ values, the speed of
sound in neutron stars at $\beta$-equilibrium almost always show a peak as
long as the $J_{0}$ characterizing the stiffness of high-density SNM EOS is
not too far above its currently known most probable value of about $-190$ MeV.
Moreover, the spinodal decomposition onset density $\rho_{t}$ strongly depends
on the high-density behavior of nuclear symmetry energy $E_{\rm sym}(\rho)$.

##  4 Summary

In summary, within a meta EOS model we studied the speed of sound and the
possible onset of high-density spinodal decomposition in both nucleonic matter
with fixed isospin asymmetries and the $n+p+e+\mu$ matter in neutron stars at
$\beta$-equilibrium. While in nucleonic matter with fixed isospin asymmetries
the speed of sound may continuously increase with density and does not undergo
a spinodal decomposition especially for neutron-rich matter if the $E_{\rm
sym}(\rho)$ is very stiff, the sound speed in neutron star matter at
$\beta$-equilibrium almost always show a peak at certain density depending
strongly on the high-density behavior of $E_{\rm sym}(\rho)$ if the skewness
$J_{0}$ of SNM EOS is not too far above its currently known most probable vale
of about $-190$ MeV. Moreover, the spinodal decomposition onset density
depends sensitively on the high-density behavior of nuclear symmetry energy
$E_{\rm sym}(\rho)$. Furthermore, in the case of using a super-soft $E_{\rm
sym}(\rho)$, the spinodal decomposition occurs at a density lower than the
critical density above which isospin separation instability occurs.

As we have stated in several places above, our work has caveats. A major one
is that all calculations are carried out with the minimum model of neutron
stars consisting of nucleons, electrons and muons only without treating the
physics of new phases and new degrees of freedom that may appear. In
particular, it is known that the appearance of Delta resonances and some
hyperons is strongly correlated with the density dependence of nuclear
symmetry energy, albeit depending also sensitively on their poorly known
couplings/interactions with other particles Cai ; Drago . Some of these
particles may start appearing at densities around $(2-3)\rho_{0}$. It thus
remains an interesting question how the appearance of these particles may
affect the density dependence of the speed of sound as well as the onset
densities or possible removals of the instabilities we studied in this work.

## Acknowledgement

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

## References

  * (1) I. Bombaci, U. Lombardo, Phys. Rev. C 44, 1892 (1991).
  * (2) B.A. Li, À. Ramos, G. Verde, I. Vidaňa (Eds.), Topical issue on nuclear symmetry energy, Eur. Phys. J. A 50, No.2 (2014).
  * (3) J.M. Lattimer, M. Prakash, Phys. Rep. 333, 121 (2000).
  * (4) S. Kubis, Phys. Rev. C 76, 025801 (2007).
  * (5) J. Xu, L.W. Chen, B.A. Li, H.R. Ma, Astrophys. J. 697, 1549 (2009).
  * (6) G.F. Bertsch, P.J. Siemens, Phys. Lett. B 126, 9 (1983).
  * (7) P.J. Siemens, Nature 305, 410 (1983).
  * (8) J.A. Lopez, P.J. Siemens, Nucl. Phys. A 431, 728 (1984).
  * (9) L. Goodman, J.I. Kapusta, A.Z. Mekjian, Phys. Rev. C 30, 851 (1984).
  * (10) B.A. Li, S. Pratt, P.J. Siemens, Phys. Rev. C 37, 728 (1988).
  * (11) H. Muller, B.D. Serot, Phys. Rev. C 52, 2072 (1995).
  * (12) L.G. Moretto, K. Tso, N. Colonna, G.J. Wozniak, Phys. Rev. Lett. 69, 1884 (1992).
  * (13) W. Bauer, G.F. Bertsch, H. Schulz, Phys. Rev. Lett. 69, 1888 (1992).
  * (14) B.A. Li, D.H.E. Gross, Nucl. Phys. A 554, 257 (1993).
  * (15) H.M. Xu, J.B. Natowitz, C.A. Gagliardi, R.E. Tribble, C.Y. Wong, W.G. Lynch, Phys. Rev. C 48, 933 (1993).
  * (16) P. Chomaz, M. Colonna, J. Randrup, Phys. Rep. 389, 263 (2004).
  * (17) M. Colonna, Prog. Part. Nucl. Phys. 113, 103775 (2020).
  * (18) P. Danielewicz, H. Lin, J.R. Stone, Y. Iwata, Phys. Lett. B 811, 135957 (2020).
  * (19) D.H. Rischke, M. Gyulassy, Nucl. Phys. A 597, 701 (1996).
  * (20) P. Danielewicz, R.A. Lacey, P.B. Gossiaux, C. Pinkenburg, P. Chung, J.M. Alexander, R.L. McGrath, Phys. Rev. Lett. 81, 2438 (1998).
  * (21) D.H. Rischke, Y. Pürsün, J.A. Maruhn, H. Stoecker, W. Greiner, Acta Phys. Hung. A 1, 309 (1995).
  * (22) B.A. Li, C.M. Ko, Phys. Rev. C 58, R1382 (1998).
  * (23) J. Steinheimer, J. Randrup, Phys. Rev. Lett. 109, 212301 (2012).
  * (24) B.A. Li, M. Magno, Phys. Rev. C 102, 045807 (2020).
  * (25) M.G. Alford, S. Han, M. Prakash, Phys. Rev. D 88, 083013 (2013).
  * (26) S. Han, AIP Conference Proceedings 2127, 020006 (2019).
  * (27) E. Annala, T. Gorda, A. Kurkela, J. Nättilä, A. Vuorinen, Nature Phys. 16, 907 (2020).
  * (28) G. Grams, R. Somasundaram, J. Margueron, E. Khan, Phys. Rev. C 106, 044305 (2022).
  * (29) H. K. Lee, Y.L. Ma, W.G. Paeng, M. Rho, Mod. Phys. Lett. A 37, 2230003 (2022).
  * (30) T. Malik, M. Ferreira, B.K. Agrawal, C. Providência, Astrophys. J. 930, 17 (2022).
  * (31) N. Lozano, V. Tran, P. Jaikumar, Galaxies 10, 79 (2022).
  * (32) L. McLerran, S. Reddy, Phys. Rev. Lett. 122, 122701 (2019).
  * (33) H. Tan, V. Dexheimer, J. Noronha-Hostler, N. Yunes, Phys. Rev. Lett. 128, 161101 (2022).
  * (34) C. Drischler, S. Han, S. Reddy, Phys. Rev. C 105, 035808 (2022).
  * (35) R.M. Aguirre, Phys. Rev. D 105, 116023 (2022).
  * (36) S. Altiparmak, C. Ecker, L. Rezzolla, arXiv:2203.14974.
  * (37) J. Braun, A. Geißel, B. Schallmo, arXiv:2206.06328.
  * (38) A. Kanakis-Pegios, P.S. Koliogiannis, C.C. Moustakidis, Phys. Rev. C 102, 055801 (2020).
  * (39) T. Kojo, AAPPS Bull. 31, 11 (2021).
  * (40) T.F. Motta, P.A.M. Guichon, A.W. Thomas, Nucl. Phys. A 1009, 122157 (2021).
  * (41) M. Marczenko, L. McLerran, K. Redlich, C. Sasaki, arXiv:2207.13059.
  * (42) D. Kumar, H. Mishra, T. Malik, arXiv:2110.00324.
  * (43) M.Z. Han, Y.J. Huang, S.P. Tang, Y.Z. Fan, arXiv:2207.13613.
  * (44) C. Ecker, L. Rezzolla, Astrophys. J. Lett. 939, L35 (2022).
  * (45) Y. Fujimoto, K. Fukushima, L.D. McLerran, M. Praszalowicz, Phys. Rev. Lett. 129, 252702 (2022).
  * (46) I. Tews, J. Carlson, S. Gandolfi, S. Reddy, Astrophys. J. 860 149 (2018).
  * (47) H. Tan, T. Dore, V. Dexheimer, J. Noronha-Hostlerm, N. Yunes, Phys. Rev. D 105, 023018 (2022).
  * (48) W.J. Xie, B.A. Li, Astrophys. J. 883, 174 (2019).
  * (49) W.J. Xie, B.A. Li, Astrophys. J. 899, 4 (2020).
  * (50) B.A. Li, L.W. Chen, C.M. Ko, Phys. Rep. 464, 113 (2008).
  * (51) C. Providência, S.S. Avancini, R. Cavagnoli, S. Chiacchiera, C. Ducoin, F. Grill, J. Margueron, D.P. Menezes, A. Rabhi, I. Vidaña, Euro. Phys. J. A 50, 44 (2014).
  * (52) B.A. Li, P.G. Krastev, D.H. Wen, N.B. Zhang, Euro. Phys. J. A 55, 117 (2019).
  * (53) M. Ferreira, A. Rabhi, C. Providência, Eur. Phys. J. A 57, 263 (2021).
  * (54) M. Kutschera, W. Wojcik, Phys. Lett. B 325, 271 (1994).
  * (55) J. Margueron, J. Navarro, N. van Giai, Phys. Rev. C 66, 014303 (2002).
  * (56) J. Rikovska Stone, J.C. Miller, R. Koncewicz, P.D. Stevenson, M.R. Strayer, Phys. Rev. C 68, 034324 (2003).
  * (57) A.Lavagno, D. Pigato, Physica A 531, 121595 (2019).
  * (58) A.R. Raduta, Phys. Lett. B 814, 136070 (2021).
  * (59) W.Z. Jiang, B.A. Li, L.W. Chen, Astrophys. J. 756, 56 (2012).
  * (60) B.J. Cai, F.J. Fattoyev, B.A. Li, W.G. Newton, Phys. Rev. C 92, 015802 (2015).
  * (61) Z. Miao, A. Li, Z. Zhu, S. Han, Astrophys. J. 904, 103 (2020).
  * (62) S.P. Tang, J.L. Jiang, W.H. Gao, Y.Z. Fan, D.M. Wei, Phys. Rev. D 103, 063026 (2021).
  * (63) A. Li, Z. Miao, S. Han, B. Zhang, Astrophys. J. 913, 27 (2021).
  * (64) R. Somasundaram, J. Margueron, Europhysics Letters 138, 14002 (2022).
  * (65) W.J. Xie, B.A. Li, Phys. Rev. C 103, 035802 (2021).
  * (66) H. D. Thi, C. Mondal and F. Gulminelli, [arXiv:2302.08482 [nucl-th]].
  * (67) J. Margueron, R. Hoffmann Casali, F. Gulminelli, Phys. Rev. C 97, 025805 (2018).
  * (68) J. Margueron, R.H. Casali, F. Gulminelli, Phys. Rev. C 97, 025806 (2018).
  * (69) C. Mondal, F. Gulminelli, Phys. Rev. C 107, no.1, 015801 (2023).
  * (70) N.B. Zhang, B.A. Li, J. Xu, Astrophysical J. 859, 90 (2018).
  * (71) N.B. Zhang, B.A. Li, Euro. Phys. J. A 55, 39 (2019).
  * (72) U. Garg, G. Colò, Prog. Part. Nucl. Phys. 101, 55 (2018).
  * (73) B.A. Li, X. Han, Phys. Lett. B  727, 276 (2013).
  * (74) M. Oertel, M. Hempel, T. Klähn, S. Typel, Rev. Mod. Phys. 89, 015007 (2017).
  * (75) B.A. Li, B.J. Cai, W.J. Xie, N.B. Zhang, Universe 7, 182 (2021).
  * (76) D. Adhikari et al. (PREX Collaboration), Phys. Rev. Lett. 126, 172502 (2021).
  * (77) B.T. Reed, F.J. Fattoyev, C.J. Horowitz, J. Piekarewicz, 126, 172503 (2021).
  * (78) P.G. Reinhard, X. Roca-Maza, W. Nazarewicz, Phys. Rev. Lett. 127, 232501 (2021).
  * (79) R. Essick, I. Tews, P. Landry, A. Schwenk, Phys. Rev. Lett. 127, 192701 (2021).
  * (80) T.G. Yue, L.W. Chen, Z. Zhang, Y. Zhou, Phys. Rev. Research 4, L022054 (2022).
  * (81) J.M. Lattimer, Particles 6, 30 (2023).
  * (82) J. Estee et al. (S$\pi$RIT Collaboration), Phys. Rev. Lett. 126, 162701 (2021).
  * (83) H.J. Xu, Acta Phys. Polon. Supp. 16, 1-A30 (2023).
  * (84) Skyy V. Pineda et al., Phys. Rev. Lett. 127, 182503 (2021).
  * (85) C. Drischler, R.J. Furnstahl, J.A. Melendez, D.R. Phillips, Phys. Rev. Lett. 125, 202702 (2020).
  * (86) D. Adhikari et al, Phys. Rev. Lett. 129, 042501 (2022).
  * (87) P.G. Reinhard, X. Roca-Maza, W. Nazarewicz, Phys. Rev. Lett. 129, 232501 (2022).
  * (88) Z. Zhang, L.W. Chen, arXiv:2207.03328 [nucl-th].
  * (89) S.Tagami, T.Wakasa, M. Yahiro, arXiv:2210.03313 [nucl-th].
  * (90) E. Yüksel and N. Paar, Phys. Lett. B 836, 137622 (2023)
  * (91) B. Biswas, Astrophys. J. 921, no.1, 63 (2021).
  * (92) M. A. Corona, M. Cadeddu, N. Cargioli, P. Finelli and M. Vorabbi, Phys. Rev. C 105, no.5, 055503 (2022)
  * (93) C. Mondal, B.K. Agrawal, J.N. De, S.K. Samaddar, M. Centelles, X. Viñas, Phys. Rev. C 96, 021302 (2017).
  * (94) M. Dutra, O. Louren, J.S.S. Martins, A. Delfino, J.R. Stone, P.D. Stevenson, Phys. Rev. C 85, 035201 (2012).
  * (95) M. Dutra, O. Lourenço, S.S. Avancini, B.V. Carlson, A. Delfino, D.P. Menezes, C. Providência, S. Typel, J.R. Stone, Phys. Rev. C 90, 055203 (2014).
  * (96) I. Tews, J.M. Lattimer, A. Ohnishi, E.E. Kolomeitsev, Astrophys. J. 848, 105 (2017).
  * (97) N.B. Zhang, B.J. Cai, B.A. Li, W.G. Newton, J. Xu, Nucl. Sci. Tech. 28, 181 (2017).
  * (98) N.B. Zhang, B.A Li, Astrophys. J. 921, 111 (2021).
  * (99) L.W. Chen, Eur. Phys. J. Web Conf. 88, 00017 (2018).
  * (100) B.J. Cai, B.A. Li, Annals of Physics 444, 169062 (2022).
  * (101) M. Kutschera, Phys. Lett. B 340, 1 (1994).
  * (102) B.A. Li, Phys. Rev. Lett. 88, 192701 (2002).
  * (103) M. Kutschera, W. Wojcik, Phys. Rev. C 47, 1077 (1993).
  * (104) A. Szmagliński, W. Wójcik, M. Kutschera., Acta Phys. Pol. B 37, 277 (2006).
  * (105) D.H. Wen, B.A. Li, L.W. Chen, Phys. Rev. Lett. 103, 211102 (2009).
  * (106) D.T. Khoa, W. von Oertzen, A.A. Ogloblin, Nucl. Phys. A 602, 98 (1996).
  * (107) D.N. Basu, P. Roy Chowdhury, C. Samanta, Nucl. Phys. A 811, 140 (2008).
  * (108) S. Banik, D. Bandyopadhyay, J. Phys. G 26, 1495 (2000).
  * (109) S. Kubis, M. Kutschera, Nucl. Phys. A 720, 189 (2003).
  * (110) A. Sorensen et al., arXiv:2301.13253 [nucl-th].
  * (111) W.J. Xie, B.A. Li, J. Phys. G 48, 025110 (2021).
  * (112) N.B. Zhang, B.A. Li, Astrophys. J. 879, 99 (2019).
  * (113) C. Ducoin, J. Margueron, C. Providencia, I. Vidana, Phys. Rev. C 83, 045810 (2011).
  * (114) R.B. Wiringa, V. Fiks, A. Fabrocini, Phys. Rev. C 38, 1010 (1988).
  * (115) J.M. Lattimer, C.J. Pethick, M. Prakash, P. Haensel, Phys. Rev. Lett. 66, 2701 (1991).
  * (116) J.M. Lattimer, M. Prakash, Astrophys. J. 550, 426 (2001).
  * (117) J.P. Blaizot, Phys. Rep. 64, 171 (1980).
  * (118) D. Lai, Mon. Not. R. Astron. Soc. 270, 611 (1994).
  * (119) A. Reisenegger, P. Goldreich, Astrophys. J. 395, 240 (1992).
  * (120) B.A. Li, C.M. Ko, Nucl. Phys. A 618, 498 (1997).
  * (121) V. Baran, A. Larionov, M. Colonna, M. Di Toro, Nucl. Phys. A 632, 287 (1998).
  * (122) L. Shi, P. Danielewicz, Euro. Phys. Lett. 49, 34 (2000).
  * (123) H.S. Xu, M.B. Tsang, T.X. Liu et al., Phys. Rev. Lett. 85, 716 (2000).
  * (124) M. Ferreira, R. Câmara Pereira, C. Providência, Phys. Rev. D 101, 123030 (2020).
  * (125) A. Drago, A. Lavagno, G. Pagliara, D. Pigato, Phys. Rev. C 90, 065809 (2014).