#  Exploring the first-order phase transition in neutron stars using the
parity doublet model and NJL-type quark model

Bikai Gao<EMAIL_ADDRESS>Department of Physics, Nagoya
University, Nagoya 464-8602, Japan  School of Frontier Sciences, Nanjing
University, Suzhou, 215163, China  Wen-Li Yuan<EMAIL_ADDRESS>School of
Physics and State Key Laboratory of Nuclear Physics and Technology, Peking
University, Beijing 100871, China  Masayasu Harada
<EMAIL_ADDRESS>Department of Physics, Nagoya University,
Nagoya 464-8602, Japan  Kobayashi-Maskawa Institute for the Origin of
Particles and the Universe, Nagoya University, Nagoya, 464-8602, Japan  Yong-
Liang Ma<EMAIL_ADDRESS>School of Frontier Sciences, Nanjing University,
Suzhou, 215163, China

###### Abstract

We investigate the possibility and impacts of a first-order phase transition
from hadronic matter to quark matter in neutron stars (NSs) using two specific
models: the parity doublet model (PDM) for the hadronic phase and the Nambu-
Jona-Lasinio (NJL) type model for the quark phase. By combining these models,
we construct hybrid equations of state (EOSs) that capture the transition
between the two phases. We explore the parameter space of both models to
identify the conditions under which a first-order phase transition can occur
and study its effects on NS properties. We identify the suitable parameter
space and constrain the onset density of the first-order phase transition. For
$m_{0}$ = 500 MeV—the chiral invariant mass in PDM, the phase transition
occurs between 1.9$n_{0}$ and 2.95$n_{0}$ and ends between 2.1$n_{0}$ and
3.6$n_{0}$. Increasing $m_{0}$ to 600 MeV shifts the phase transition to
higher densities, occurring between 2.9$n_{0}$ and 4.1$n_{0}$ and ending
between 3.4$n_{0}$ and 4.6$n_{0}$.

##  I Introduction

Quantum Chromodynamics (QCD) is believed to be the fundamental theory that
describes the strong interactions between quarks and gluons. Under normal
conditions, quarks and gluons are confined within hadrons, such as protons,
neutrons, and pions. These composite particles serve as effective degrees of
freedom (d.o.f.) for describing the physics in vacuum and low
density/temperature. However, the situation changes with the increasing of
temperature or the density as hadrons begin to overlap and quarks and gluons
d.o.f. will come into play. Such extreme forms of matter are believed to have
existed in the early universe and are present in astrophysical objects like
neutron stars (NSs) (see, e.g., Refs. Fukushima and Hatsuda (2011); Fukushima
and Sasaki (2013); Baym  _et al._ (2018); Ma and Rho (2020); Annala  _et al._
(2020); Lovato  _et al._ (2022); Brandes and Weise (2024) and references
therein).

Along the temperature axis in hot QCD, a heated hadron resonance gas smoothly
transforms into a quark-gluon plasma, despite the apparent differences between
these two phases. Experiments have successfully created such hot matter
through high-energy heavy ion collisions, where the colliding energy is
converted into heat. A coherent picture of hot QCD matter has emerged through
the combined efforts of experimental analyses, ab-initio lattice simulations
of QCD, and model calculations Fukushima and Hatsuda (2011); Lovato  _et al._
(2022); Fukushima (2012); Petreczky (2012).

However, the situation is much different for cold dense matter, especially
about the exact nature of the transition from nuclear matter to quark matter.
Theoretically, this cannot be accessed due to the lack of solid approach
anchored on nonperturbative QCD. Moreover, the model dependent predictions
cannot be diagnosed by terrestrial experiments since such high density cannot
be produced or Lattice QCD simulation due to the sign problem. The only
constraints on nuclear matter properties are from the experimental data around
the saturation density $n_{0}$ ($n_{0}=0.16$ fm3) and predictions from the
perturbative QCD at extreme high density $\mathrel{\hbox to0.0pt{\lower
3.0pt\hbox{\hskip 1.0pt$\sim$}\hss}\raise 1.0pt\hbox{$>$}}50n_{0}$.

Since the detection of the gravitational waves (GWs) from a binary NS merger
GW170817 Abbott (2017); Abbott  _et al._ (2017, 2018), astrophysics enters the
multi-messenger era. In combination with the observation of massive pulsars
around two-solar mass Miller  _et al._ (2021); Riley  _et al._ (2021);
Vinciguerra  _et al._ (2024), NS has become a natural laboratory for
investigating the properties of matter under extreme conditions, particularly
at ultra-high densities up to $\sim 10n_{0}$. The characteristics of neutron
stars are determined by the equation of state (EOS) of the matter composing
them. Therefore, the theoretically obtained EOSs can be constrained by
observational data, such as the mass-radius ($M\text{-}R$) relation of a NS,
the tidal deformability, and the signals of gravitational waves stemming from
binary NS mergers Annala  _et al._ (2020); Brandes and Weise (2024);
Hotokezaka  _et al._ (2011, 2013); De Pietri  _et al._ (2016); Maione  _et
al._ (2016); Moustakidis  _et al._ (2017); Radice  _et al._ (2017);
Paschalidis  _et al._ (2018); Annala  _et al._ (2018); Torres-Rivas  _et al._
(2019); Most  _et al._ (2019); Tews  _et al._ (2018); Tsang  _et al._ (2018);
Bauswein  _et al._ (2019); Orsaria  _et al._ (2019); Most  _et al._ (2020);
Yang  _et al._ (2021); Sekiguchi  _et al._ (2011); Choi  _et al._ (2021);
Annala  _et al._ (2022); Raaijmakers  _et al._ (2021); Prakash  _et al._
(2021, 2021); Kedia  _et al._ (2022); Fujimoto  _et al._ (2023); Guo  _et al._
(2023).

One of the key questions in the study of NS EOSs is whether there is a
possible transition from hadronic matter to quark matter—either a phase
transition or a crossover. Such a transition could have significant
implications for the structure and properties of NSs (see, e.g., Baym  _et
al._ (2018); Lenzi and Lugones (2012); Kojo (2021); Benic  _et al._ (2015)).
When the transition is a first-order phase transition, there can be a sudden
release or absorption of energy, leading to changes in the equilibrium state
of the star. This can manifest as a sudden change in the radius of the neutron
star, as the internal pressure and composition adjust to the new equilibrium
state.

In previous works, by assuming the transition is a cossover transition, some
of us constructed unified EOSs by interpolating the hadronic EOS from a
hadronic model based on the parity doublet structure and the quark EOS from an
Nambu-Jona-Lasinio (NJL)-type quark model Minamikawa  _et al._ (2021a, b); Gao
_et al._ (2022); Kong  _et al._ (2023); Gao  _et al._ (2024a). In this paper,
we devote ourselves to investigate the possibility and the impacts of the
first-order phase transition in dense matter using the parity doublet model
(PDM) for the hadronic phase and the NJL-type model for the quark phase.

The PDM offers a unique perspective of baryon structure by incorporating
chiral symmetry and its restoration, allowing the parity partners to have a
degenerate chiral invariant mass $m_{0}$ at high densities when chiral
symmetry is restored. The notion of chiral invariant mass was first introduced
in DeTar and Kunihiro (1989); Jido  _et al._ (2001) and is consistent with the
recent lattice QCD simulation Aarts  _et al._ (2015, 2017, 2019) as well as
the skyrmion crystal approach to nuclear physics Ma  _et al._ (2013, 2014).
Effective models based on parity doublet structure is widely used for studying
the hadron structure and nuclear/NS physicsNishihara and Harada (2015);
Minamikawa  _et al._ (2023a); Hatsuda and Prakash (1989); Zschiesche  _et al._
(2007); Dexheimer  _et al._ (2008); Sasaki and Mishustin (2010); Sasaki  _et
al._ (2011); Gallas  _et al._ (2011); Steinheimer  _et al._ (2011); Paeng  _et
al._ (2012); Dexheimer  _et al._ (2013); Paeng  _et al._ (2013); Mukherjee
_et al._ (2017a); Suenaga (2018); Takeda  _et al._ (2018); Abuki  _et al._
(2018); Marczenko  _et al._ (2019a); Motohiro  _et al._ (2015); Yamazaki and
Harada (2019); Mukherjee  _et al._ (2017b); Gao  _et al._ (2022); Marczenko
_et al._ (2022); Minamikawa  _et al._ (2023b); Baym  _et al._ (2018);
Minamikawa  _et al._ (2021a, b); Marczenko  _et al._ (2019b); Kong  _et al._
(2023); Gao  _et al._ (2024b) . The key feature of the PDM is that a larger
$m_{0}$ leads to a weaker $\sigma$ coupling to nucleon since nucleon does not
have to acquire its mass entirely from the $\sigma$ field. The corresponding
$\omega$ field also becomes smaller due to the equilibrium state at the
saturation density $n_{0}$. Then at density larger than $n_{0}$, the $\sigma$
field decreases while the $\omega$ field increases. As a result, a larger
$m_{0}$ would lead to a softer EOS.

By combining the PDM and the NJL-type model, we aim to construct a hybrid EOS
that captures the transition from hadronic matter to quark matter. We will
explore the parameter space of both models to identify the conditions under
which a first-order phase transition can occur. The resulting EOSs will be
compared with observational constraints from NS observation to assess their
compatibility with the latest astronomical data.

This paper is organized as follows. In Sec. II, we present the EOSs used in
our study. We describe the hadronic matter EOS based on the PDM and the quark
matter EOS derived from the NJL-type model. We also discuss the construction
of the hybrid EOS incorporating a first-order phase transition. In Sec.III, we
investigate the neutron star properties using the EOSs constructed in Sec.II,
focusing on the mass-radius relation and the constraints on the onset density
of the phase transition. Finally, we summarize our findings and discuss their
implications in Sec.IV.

##  II EQUATION OF STATE

In this section, we develop EOSs for neutron star matter in the hadron and
quark phases separately, and then explore the unified EOS with a first-order
phase transition by integrating them.

###  II.1 NUCLEAR MATTER EOS

Through out this work, we assume that the hyperons do not enter into the
matter in hadronic phase and consider the case $N_{f}=2$. Following Ref.
Motohiro  _et al._ (2015), we express the thermodynamic potential based on
parity doublet structure as

$\displaystyle\Omega_{\mathrm{PDM}}=$ $\displaystyle
V\left(\sigma\right)-V\left(\sigma_{0}\right)-\frac{1}{2}m_{\omega}^{2}\omega^{2}-\frac{1}{2}m_{\rho}^{2}\rho^{2}$
(1)
$\displaystyle-\lambda_{\omega\rho}\left(g_{\omega}\omega\right)^{2}\left(g_{\rho}\rho\right)^{2}$
$\displaystyle-2\sum_{i=+,-}\sum_{\alpha=p,n}\int^{k_{f}}\frac{\mathrm{d}^{3}\mathbf{p}}{(2\pi)^{3}}\left(\mu_{\alpha}^{*}-E_{\mathrm{p}}^{i}\right).$

Here $i=+,-$ denote for the parity of nucleon—since in PDM, not only the
positive parity baryons, proton and neutron, but also the negative parity
baryons which are regarded as the parity partner of positive parity baryons
are involved—and $E_{{\bf p}}^{i}=\sqrt{{\bf p}^{2}+m_{i}^{2}}$ is the energy
of baryon with mass $m_{i}$ and momentum ${\bf p}$. The potential $V(\sigma)$
is taken as

$\displaystyle
V(\sigma)=-\frac{1}{2}\bar{\mu}^{2}\sigma^{2}+\frac{1}{4}\lambda_{4}\sigma^{4}-\frac{1}{6}\lambda_{6}\sigma^{6}-m_{\pi}^{2}f_{\pi}\sigma.$
(2)

The total thermodynamic potential for the matter is obtained by including the
effects of leptons as

$\displaystyle\Omega_{{\rm H}}=\Omega_{{\rm PDM}}+\sum_{l=e,\mu}\Omega_{l},$
(3)

where $\Omega_{l}(l=e,\mu)$ are the thermodynamic potentials for leptons,

$\Omega_{l}=-2\int^{k_{F}}\frac{d^{3}\mathbf{p}}{(2\pi)^{3}}\left(\mu_{l}-E_{\mathbf{p}}^{l}\right).$
(4)

The parameters in the PDM are determined by fitting the nuclear matter
properties around saturation density $n_{0}$ and the pion decay constant for
different $m_{0}$. Using the explicit values evaluated in Ref. Gao  _et al._
(2024a) we can then calculate the corresponding EOS shown in Fig. 1 in the
hadronic phase with different choices of chiral invariant mass $m_{0}$.

Figure 1: EOS for different values of the chiral invariant mass $m_{0}$.

From Fig. 1, one easily identify that a larger chiral invariant mass $m_{0}$—a
feature of the PDM—leads to a relatively softer EOS. In addition, one can also
observe a transition due to the entering of the nucleon with negative
parity—another feature of the PDM, as shown in Fig. 2. In a subsequent
section, we will employ the Maxwell construction to derive the EOS with a
first-order phase transition within the PDM, which is associated with the
inclusion of the negative parity nucleon.

Figure 2: Number density for different components with different values of the
chiral invariant mass $m_{0}$.

The onset density of the transition, characterized by the appearance of the
negative parity nucleon N(1535), is strongly influenced by the choice of the
chiral invariant mass $m_{0}$. In the PDM, the nucleon N(939) and its chiral
partner N(1535) acquire their masses through two different mechanisms: the
chiral invariant mass $m_{0}$, which is independent of the chiral condensate,
and the mass generated by the chiral condensate, which is proportional to the
sigma field $\sigma$. The mass difference between N(939) and N(1535) can be
expressed as:

$\displaystyle\Delta m=g_{1}\sigma+g_{2}\sigma=(g_{1}+g_{2})\sigma$ (5)

where $g_{1}$ and $g_{2}$ are coupling constants. As the matter density
increases, the chiral condensate $\sigma$ decreases due to the partial
restoration of chiral symmetry. For a smaller value of $m_{0}$, the coupling
constants $g_{1}$ and $g_{2}$ must be larger to reproduce the physical masses
of N(939) and N(1535) in vacuum111The coupling constants $g_{1,2}$ are
calculated as  $\displaystyle
g_{1,2}=\frac{1}{2f_{\pi}}\left(\sqrt{(m_{-}+m_{+})^{2}-4m_{0}^{2}}\pm(m_{-}-m_{+})\right).$
We can easily observe that smaller $m_{0}$ lead to larger value of $g_{1,2}$..
In that case, the mass difference $\Delta m$ is more sensitive to changes in
the chiral condensate $\sigma$ when $m_{0}$ is small since the coupling
constant $g_{1}$ and $g_{2}$ are comparatively larger. As a result, for a
smaller $m_{0}$, the mass of N(1535) approaches the mass of N(939) more
rapidly as the density increases. This means that N(1535) becomes
energetically favorable and appears in the system at a lower density compared
to the case with a larger $m_{0}$. This is confirmed in the explicit
calculation shown in Fig. 2 that a smaller chiral invariant mass $m_{0}$ leads
to an earlier emergence of the negative parity baryon N(1535).

###  II.2 QUARK MATTER EOS

Following Refs.Baym  _et al._ (2018, 2019), we use the NJL-type model to
describe the quark matter. The model includes three-flavors and U(1)A anomaly
effects through the quark version of the Kabayashi-Maskawa-’t Hooft (KMT)
interaction Kobayashi and Maskawa (1970); ’t Hooft (1986). The coupling
constants are chosen to be the Hatsuda-Kunihiro parameters which successfully
reproduce the hadron phenomenology at low energy Baym  _et al._ (2018);
Hatsuda and Kunihiro (1994): $G\Lambda^{2}=1.835,K\Lambda^{5}=9.29$ with
$\Lambda=631.4\,\rm{MeV}$, see the definition below. The couplings $g_{V}$ and
$H$ characterize the strength of the vector repulsion and attractive diquark
correlations whose range will be examined later when we discuss the NS
constraints.

We can then write down the thermodynamic potential as

$\displaystyle\Omega_{\mathrm{CSC}}=$
$\displaystyle\,\Omega_{s}-\Omega_{s}\left[\sigma_{f}=\sigma_{f}^{0},d_{j}=0,\mu_{q}=0\right]$
(6)
$\displaystyle+\Omega_{c}-\Omega_{c}\left[\sigma_{f}=\sigma_{f}^{0},d_{j}=0\right],$

where the subscript “0” stands for the values in vacuum, and

$\displaystyle\Omega_{s}=-2\sum_{i=1}^{18}\int^{\Lambda}\frac{d^{3}\mathbf{p}}{(2\pi)^{3}}\frac{\epsilon_{i}}{2},$
(7)
$\displaystyle\Omega_{c}=\sum_{i=u,d,s}\left(2G\sigma_{i}^{2}+Hd_{i}^{2}\right)-4K\sigma_{u}\sigma_{d}\sigma_{s}-g_{V}n_{q}^{2},$
(8)

with $\sigma_{f}$ being the chiral condensates of quark flavor $f$, $d_{j}$
being quark condensates, and $n_{q}$ being the quark density. In Eq.(7),
$\epsilon_{i}$ are energy eigenvalues obtained from inverse propagator in
Nambu-Gorkov bases

$S^{-1}(k)=\left(\begin{array}[]{lc}\gamma_{\mu}k^{\mu}-\hat{M}+\gamma^{0}\hat{\mu}&\gamma_{5}\sum_{i}\Delta_{i}R_{i}\\\
-\gamma_{5}\sum_{i}\Delta_{i}^{*}R_{i}&\gamma_{\mu}k^{\mu}-\hat{M}-\gamma^{0}\hat{\mu}\end{array}\right),$
(9)

where

$\displaystyle
M_{i}=m_{i}-4G\sigma_{i}+K\left|\epsilon_{ijk}\right|\sigma_{j}\sigma_{k},$
(10) $\displaystyle\Delta_{i}=-2Hd_{i},$
$\displaystyle\hat{\mu}=\mu_{q}-2g_{V}n_{q}+\mu_{3}\lambda_{3}+\mu_{8}\lambda_{8}+\mu_{Q}Q,$
$\displaystyle(R_{1},R_{2},R_{3})=(\tau_{7}\lambda_{7},\tau_{5}\lambda_{5},\tau_{2}\lambda_{2}).$

$S^{-1}(k)$ is $72\times 72$ matrix in terms of the color, flavor, spin, and
Nambu-Gorkov basis, which has 72 eigenvalues. $M_{u,d,s}$ are the constituent
masses of $u,d,s$ quarks and $\Delta_{1,2,3}$ are the gap energies. The
$\mu_{3,8}$ are the color chemical potentials which will be tuned to achieve
the color neutrality. The total thermodynamic potential including the effect
of leptons is

$\Omega_{\mathrm{Q}}=\Omega_{\mathrm{CSC}}+\sum_{l=e,\mu}\Omega_{l}.$ (11)

The mean fields are determined from the gap equations,

$0=\frac{\partial\Omega_{\mathrm{Q}}}{\partial\sigma_{i}}=\frac{\partial\Omega_{\mathrm{Q}}}{\partial
d_{i}}.$ (12)

From the conditions for electromagnetic charge neutrality and color charge
neutrality, we have

$n_{j}=-\frac{\partial\Omega_{\mathrm{Q}}}{\partial\mu_{j}}=0,$ (13)

where $j=3,8,Q$. The baryon number density $n_{B}$ is determined as

$n_{q}=-\frac{\partial\Omega_{\mathrm{Q}}}{\partial\mu_{q}},$ (14)

where $\mu_{q}$ is $1/3$ of the baryon number chemical potential. After
determined all the values, we obtain the pressure as

$P_{\mathrm{Q}}=-\Omega_{\mathrm{Q}}.$ (15)

###  II.3 FIRST-ORDER PHASE TRANSITION

In this section, we construct the EOS for NS matter incorporating a first-
order phase transitionLenzi and Lugones (2012); Benic  _et al._ (2015). The
hadronic phase is described by the PDM, which includes the parameter $m_{0}$,
representing the chiral invariant mass. The quark phase is modeled using an
NJL-type quark model with two key parameters: $H$, characterizing the strength
of attractive diquark correlations, and $g_{V}$, denoting the strength of
vector repulsion. By adjusting the values of these parameters, we can control
the stiffness of the quark matter EOS and, consequently, the onset and
properties of the phase transition.

To construct a hybrid EOS with a first-order phase transition, we first fix
the value of $m_{0}$ in the PDM. As an example, we set $m_{0}=500$ MeV, which
corresponds to a relatively stiff hadronic EOS. We then explore various
combinations of $(H,g_{V})$ in the NJL-type quark model to find an
intersection point between the hadronic and quark matter EOSs, as shown in
Fig. 3(a). The intersection point indicates the baryon chemical potential at
which the phase transition occurs, ensuring that the pressure and chemical
potential of both phases are equal, satisfying the Maxwell construction Zdunik
and Haensel (2013); Alford  _et al._ (2013) for a first-order phase
transition. The resulting hybrid EOS is presented in Fig. 3(b). If the
transition density from hadronic matter to quark matter is larger than density
where chiral phase transition—the degeneracy of the positive and negative
parity baryons—can happen, then two distinct plateaus are observed in the EOS
for $(H,g_{V})=(1.3,0.5),(1.3,0.9),(1.3,1.1)$. The first plateau, located
around $2n_{0}$, corresponds to the emergence of the negative parity baryon
$N(1535)$, a characteristic feature of the PDM. This transition is associated
with a softening of the EOS due to the partial restoration of chiral symmetry.
The second plateau represents the first-order phase transition from hadronic
matter to quark matter. The onset of this transition is determined by the
chosen values of $(H,g_{V})$ in the NJL-type quark model. During the phase
transition, the pressure remains constant while the baryon density experiences
a discontinuous jump, indicating a coexistence region between hadronic and
quark matter phases. On the other hand, if the quark-hadron transition density
is smaller than the chiral phase transition, there is only one plateau as the
blue curve $(H,g_{V})=(1.3,0.5)$ shows.

The presence of these two transitions in the hybrid EOS has significant
implications for the structure and properties of neutron stars. The softening
of the EOS due to the transitions can impact the maximum mass and radius of
neutron stars, as well as their tidal deformability and other observable
properties. By systematically varying the parameters $m_{0}$ in the PDM and
$(H,g_{V})$ in the NJL-type quark model, we can explore a range of possible
hybrid EOSs with first-order phase transition and assess their compatibility
with observational constraints from NSs.

Figure 3: (a)Pressure as a function of baryon chemical potential for fixed
value of $m_{0}$ in the PDM and several choices of ($H,g_{V}$) combination in
the NJL-type quark model (upper panel). (b)Pressure as a function of baryon
number density in unit of normal nuclear density $n_{0}$ constructed from (a).

##  III Neutron star properties

In this section, we study the NS properties by using the EOSs constructed in
the previous section. By solving the Tolman-Oppenheimer-Volkoff (TOV) equation
for spherically symmetric and static stars, we obtain the NS mass-radius (M-R)
relation.

We first show the M-R relation in pure hadronic PDM in Fig. 4. The green and
red contours represent the NS observational constraints from LIGO-VIRGO Abbott
(2017); Abbott  _et al._ (2017, 2018) and the Neutron Star Interior
Composition Explorer (NICER) Miller  _et al._ (2021); Riley  _et al._ (2021);
Vinciguerra  _et al._ (2024), respectively. We observe that for $m_{0}$ = 700
and 800 MeV, the maximum mass of the NS is smaller than 2$M_{\odot}$
constraint. Since the first-order phase transition soften the EOS and lower
the maximum mass of the NS, it is impossible to realize the NS with
$M/M_{\odot}\geq 2$ for $m_{0}\geq 700$ MeV. Consequently, we focus on the
cases with $m_{0}$ = 500 and 600 MeV in the subsequent analysis.

Figure 4: Mass-radius relation in PDM for different values of $m_{0}$.

Next, we examine the impact of the first-order phase transition on the M-R
relation by combining the PDM with the NJL-type model. We show some typical
example in Fig. 5 and Fig. 6, with the black line representing the pure PDM
and different colors corresponding to various NJL parameter combinations
$(H,g_{V})$. For the purpose to investigate the effects of parameters $H$ and
$g_{V}$, we change one values with keeping the others fixed, as shown in Fig.
5(a) and Fig. 5(b). The cross marks indicate the points of maximum mass for
each combination.

(a)

(b)

Figure 5: Mass-radius relation for $m_{0}=500$ MeV combined with different
choices of NJL parameters and corresponding mass-center density relation. In
(a), we fix the $H$ value and change the $g_{V}$ parameter while in (b) we fix
the $g_{V}$ parameter while change the $H$.

From Fig. 5, one can conclude that the increasing of $g_{V}$ or the decreasing
of $H$ increases the onset density of the first-order phase transition. When
the phase transition occurs at a very low density, as indicated by the orange
line in Fig. 5, the magnitude of chiral symmetry breaking is suppressed. As
the onset density of the first-order phase transition is pushed to higher
densities by increasing $g_{V}$ or decreasing $H$, the chiral restoration
around 2$n_{0}$ for $m_{0}=500$ MeV begins to manifest. In the figure, the
orange, red, and blue lines represent stable NS configurations, with the
maximum mass $\sim 2m_{\odot}$ appearing after the first-order phase
transition. However, the brown line and the magenta line in Fig. 5(b) exhibit
a significant reduction in radius and mass following the phase transition,
implying an NS with a unstable quark core since the ${\rm d}M/{\rm d}n_{B}<0$
after the phase transition Alford and Sedrakian (2017); Glendenning and
Kettner (2000).

Figure 6: Mass-radius relation and mass-center density relation for
$m_{0}=600$ MeV with different NJL parameters.

We also show the case for $m_{0}=600$ MeV in Fig. 6 with some typical choices
of the NJL parameters. In this case, the onset of $N(1535)$ point is shifted
to around $3n_{0}$ and a softer quark EOS is required to form a stable NS
compared with the case in $m_{0}=500$ MeV. In the two figures in Fig. 6, we
find that the onset density of the phase transition, governed by the interplay
between hadronic and quark matter parameters, significantly influences the
stability of the NS configuration. A phase transition occurring at lower
densities allows for stable NS configurations, while a phase transition at
excessively high densities can lead to an unstable NS, leads to ${\rm d}M/{\rm
d}n_{B}<0$ after the phase transition. This is because if the onset density is
very high, only NSs with very high central densities will undergo this phase
transition. These stars are typically very massive and close to the TOV limit
(the maximum mass limit for a stable neutron star). At such high densities,
the EoS change due to the phase transition can make the core more compressible
or softer, which can destabilize the star. This destabilization may lead to
the collapse of the NS star into a black hole if the pressure support is
insufficient to counteract the gravitational pull.

Finally, by fine-tuning the NJL parameters, we can identify the suitable
parameter space and constrain the onset density of the first-order phase
transition. In Fig. 7, we show all possible combinations of $(H,g_{V})$ for
$m_{0}$ = 500, 600 and 700 MeV. The color of each point indicates the maximum
mass of the NS obtained from the corresponding parameters, as shown by the
vertical bar on the right side of each figure. A cross mark indicates that the
EOS constructed from the combination of $(H,g_{V})$ does not have an
intersection point with the EOS constructed from the PDM, while a colored
indicates that the combination is allowed. For $m_{0}$=700 MeV, all the
combination of NJL parameters lead to a NS maximum mass smaller than
2$M_{\odot}$, leading to the conclusion that the $m_{0}$ should be smaller
than 700 MeV while considering the first-order phase transition.

After considering NS with a stable quark core, we further constrain the
parameter space to the region surrounded by red lines which represents the
combinations that can simultaneously satisfy the maximum mass constraints and
the stability condition. Using these parameters, we can determine the possible
range of densities at which the first-order phase transition occurs.

We choose $m_{0}$ to be 500, 550, and 600 MeV and examine the density range of
the first-order phase transition as shown in Fig. 8. The blue line in the
figure represents the density range at which the phase transition begins,
while the red line represents the density range where the phase transition
ends. As the value of $m_{0}$ increases, we find that both the onset density
and the ending density of the phase transition increase. For $m_{0}$ = 500
MeV, the phase transition occurs between $1.9n_{0}$ and $2.95n_{0}$ and ends
between $2.1n_{0}$ and $3.6n_{0}$. When $m_{0}$ = 600 MeV, the phase
transition takes place between $2.9n_{0}$ and $4.1n_{0}$ and concludes between
$3.4n_{0}$ and $4.6n_{0}$.

(a)

(b)

(c)

Figure 7: Possible $(H,g_{V})$ combination for $m_{0}=500,600,700$ MeV. Cross
mark indicates that the EOS constructed from the combination of ($H,g_{V}$) do
not have intersect point with EOS constructed from PDM. Circle indicates that
the combination is allowed. The color shows the maximum mass of NS obtained
from the corresponding parameters, as indicated by a vertical bar at the right
side of each figure. The part surrounded by red lines represents the
combination can satisfy the maximum mass constraints and the stable condition
at the same time. Figure 8: Density range of possible first-order phase
transition as a function of $m_{0}$.

##  IV Summary and Discussion

In this paper, we have investigated the possibility and impacts of a first-
order phase transition from hadronic matter to quark matter in NSs using the
PDM for the hadronic phase and the NJL-type model for the quark phase. By
constructing hybrid EOSs and analyzing the M-R relations, we have explored the
role of key parameters in both models, such as the chiral invariant mass
($m_{0}$) in the PDM and the attractive diquark correlation strength ($H$) and
vector repulsion strength ($g_{V}$) in the NJL-type quark model.

Our findings show that the onset density of the first-order phase transition
is strongly influenced by the interplay between hadronic and quark matter
parameters. A smaller chiral invariant mass $m_{0}$ in the PDM leads to an
earlier onset of the chiral phase transition associated with the emergence of
the negative parity baryon. The stability of the NS configuration is sensitive
to the onset density of the first-order phase transition, with transitions
occurring at lower densities allowing for stable configurations, while
transitions at excessively high densities can lead to NS with an unstable
quark core. By fine-tuning the parameters $(H,g_{V})$ and considering the
stability of the NS, we have identified suitable parameter spaces and
constrained the onset density of the first-order phase transition. For $m_{0}$
= 500 MeV, the phase transition occurs between 1.9$n_{0}$ and 2.95$n_{0}$ and
ends between 2.1$n_{0}$ and 3.6$n_{0}$. For $m_{0}$ = 600 MeV, the phase
transition occurs between 2.9$n_{0}$ and 4.1$n_{0}$ and ends between
3.4$n_{0}$ and 4.6$n_{0}$.

Finally, we compare cases assuming a crossover transition with those involving
a first-order phase transition, as illustrated in Fig. 9. For each $m_{0}$
value, we select NJL parameters from the region surrounded by red lines in
Fig. 7 that yield both the largest and smallest maximum masses. Solid lines
represent cases with first-order phase transitions, while dotted lines of the
same color indicate cases assuming crossover transitions. Our analysis reveals
several key differences between these two scenarios. We observe that the
crossover assumption leads to a smooth change in the $M$-$R$ curve, while the
first-order phase transition results in a more abrupt change, often
characterized by a ”kink” in the M-R relation. Furthermore, the maximum
neutron star mass in the first-order phase transition scenario is typically
smaller than in the crossover case, which can be attributed to the sudden
softening of the equation of state during a first-order transition.

Figure 9: Comparison of the 1st order phase transition with the crossover
transition. Left panel shows pressure as the function of density and right
panel shows the corresponding $M$-$R$ curve for $m_{0}=500,600$ MeV with
different choices of NJL parameters. Solid curves are for the 1st order phase
transition and the dotted curves are for the crossover transition.

## Acknowledgments

The work of B.G., and M.H. are supported in part by JSPS KAKENHI Grant Nos.
23H05439 and 24K07045. B.G. is also supported by JST SPRING, Grant No.
JPMJSP2125. B.G. would like to take this opportunity to thank the
“Interdisciplinary Frontier Next-Generation Researcher Program of the Tokai
Higher Education and Research System.” The work of Y. L. M. was supported in
part by the National Key R&D Program of China under Grant No. 2021YFC2202900
and the National Science Foundation of China (NSFC) under Grant No. 12347103
and No. 11875147.

## References

  * Fukushima and Hatsuda (2011) K. Fukushima and T. Hatsuda, Rept. Prog. Phys. 74, 014001 (2011), arXiv:1005.4814 [hep-ph] .
  * Fukushima and Sasaki (2013) K. Fukushima and C. Sasaki, Prog. Part. Nucl. Phys. 72, 99 (2013), arXiv:1301.6377 [hep-ph] .
  * Baym  _et al._ (2018) G. Baym, T. Hatsuda, T. Kojo, P. D. Powell, Y. Song, and T. Takatsuka, Rept. Prog. Phys. 81, 056902 (2018), arXiv:1707.04966 [astro-ph.HE] .
  * Ma and Rho (2020) Y.-L. Ma and M. Rho, Prog. Part. Nucl. Phys. 113, 103791 (2020), arXiv:1909.05889 [nucl-th] .
  * Annala  _et al._ (2020) E. Annala, T. Gorda, A. Kurkela, J. Nättilä, and A. Vuorinen, Nature Phys. 16, 907 (2020), arXiv:1903.09121 [astro-ph.HE] .
  * Lovato  _et al._ (2022) A. Lovato _et al._ , (2022), arXiv:2211.02224 [nucl-th] .
  * Brandes and Weise (2024) L. Brandes and W. Weise, Symmetry 16, 111 (2024), arXiv:2312.11937 [nucl-th] .
  * Fukushima (2012) K. Fukushima, J. Phys. G 39, 013101 (2012), arXiv:1108.2939 [hep-ph] .
  * Petreczky (2012) P. Petreczky, J. Phys. G 39, 093002 (2012), arXiv:1203.5320 [hep-lat] .
  * Abbott (2017) B. P. Abbott (LIGO Scientific Collaboration and Virgo Collaboration), Phys. Rev. Lett. 119, 161101 (2017).
  * Abbott  _et al._ (2017) B. P. Abbott _et al._ (LIGO Scientific, Virgo, Fermi GBM, INTEGRAL, IceCube, AstroSat Cadmium Zinc Telluride Imager Team, IPN, Insight-Hxmt, ANTARES, Swift, AGILE Team, 1M2H Team, Dark Energy Camera GW-EM, DES, DLT40, GRAWITA, Fermi-LAT, ATCA, ASKAP, Las Cumbres Observatory Group, OzGrav, DWF (Deeper Wider Faster Program), AST3, CAASTRO, VINROUGE, MASTER, J-GEM, GROWTH, JAGWAR, CaltechNRAO, TTU-NRAO, NuSTAR, Pan-STARRS, MAXI Team, TZAC Consortium, KU, Nordic Optical Telescope, ePESSTO, GROND, Texas Tech University, SALT Group, TOROS, BOOTES, MWA, CALET, IKI-GW Follow-up, H.E.S.S., LOFAR, LWA, HAWC, Pierre Auger, ALMA, Euro VLBI Team, Pi of Sky, Chandra Team at McGill University, DFN, ATLAS Telescopes, High Time Resolution Universe Survey, RIMAS, RATIR, SKA South Africa/MeerKAT), Astrophys. J. Lett. 848, L12 (2017), arXiv:1710.05833 [astro-ph.HE] .
  * Abbott  _et al._ (2018) B. P. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. Lett. 121, 161101 (2018), arXiv:1805.11581 [gr-qc] .
  * Miller  _et al._ (2021) M. C. Miller _et al._ , Astrophys. J. Lett. 918, L28 (2021), arXiv:2105.06979 [astro-ph.HE] .
  * Riley  _et al._ (2021) T. E. Riley _et al._ , Astrophys. J. Lett. 918, L27 (2021), arXiv:2105.06980 [astro-ph.HE] .
  * Vinciguerra  _et al._ (2024) S. Vinciguerra _et al._ , Astrophys. J. 961, 62 (2024), arXiv:2308.09469 [astro-ph.HE] .
  * Hotokezaka  _et al._ (2011) K. Hotokezaka, K. Kyutoku, H. Okawa, M. Shibata, and K. Kiuchi, Phys. Rev. D 83, 124008 (2011), arXiv:1105.4370 [astro-ph.HE] .
  * Hotokezaka  _et al._ (2013) K. Hotokezaka, K. Kiuchi, K. Kyutoku, T. Muranushi, Y.-i. Sekiguchi, M. Shibata, and K. Taniguchi, Phys. Rev. D 88, 044026 (2013), arXiv:1307.5888 [astro-ph.HE] .
  * De Pietri  _et al._ (2016) R. De Pietri, A. Feo, F. Maione, and F. Löffler, Phys. Rev. D 93, 064047 (2016), arXiv:1509.08804 [gr-qc] .
  * Maione  _et al._ (2016) F. Maione, R. De Pietri, A. Feo, and F. Löffler, Class. Quant. Grav. 33, 175009 (2016), arXiv:1605.03424 [gr-qc] .
  * Moustakidis  _et al._ (2017) C. C. Moustakidis, T. Gaitanos, C. Margaritis, and G. A. Lalazissis, Phys. Rev. C 95, 045801 (2017), [Erratum: Phys.Rev.C 95, 059904 (2017)], arXiv:1608.00344 [nucl-th] .
  * Radice  _et al._ (2017) D. Radice, S. Bernuzzi, W. Del Pozzo, L. F. Roberts, and C. D. Ott, Astrophys. J. Lett. 842, L10 (2017), arXiv:1612.06429 [astro-ph.HE] .
  * Paschalidis  _et al._ (2018) V. Paschalidis, K. Yagi, D. Alvarez-Castillo, D. B. Blaschke, and A. Sedrakian, Phys. Rev. D 97, 084038 (2018), arXiv:1712.00451 [astro-ph.HE] .
  * Annala  _et al._ (2018) E. Annala, T. Gorda, A. Kurkela, and A. Vuorinen, Phys. Rev. Lett. 120, 172703 (2018), arXiv:1711.02644 [astro-ph.HE] .
  * Torres-Rivas  _et al._ (2019) A. Torres-Rivas, K. Chatziioannou, A. Bauswein, and J. A. Clark, Phys. Rev. D 99, 044014 (2019), arXiv:1811.08931 [gr-qc] .
  * Most  _et al._ (2019) E. R. Most, L. J. Papenfort, V. Dexheimer, M. Hanauske, S. Schramm, H. Stöcker, and L. Rezzolla, Phys. Rev. Lett. 122, 061101 (2019), arXiv:1807.03684 [astro-ph.HE] .
  * Tews  _et al._ (2018) I. Tews, J. Margueron, and S. Reddy, Phys. Rev. C 98, 045804 (2018), arXiv:1804.02783 [nucl-th] .
  * Tsang  _et al._ (2018) C. Y. Tsang, M. B. Tsang, P. Danielewicz, W. G. Lynch, and F. J. Fattoyev, (2018), arXiv:1807.06571 [nucl-ex] .
  * Bauswein  _et al._ (2019) A. Bauswein, N.-U. F. Bastian, D. B. Blaschke, K. Chatziioannou, J. A. Clark, T. Fischer, and M. Oertel, Phys. Rev. Lett. 122, 061102 (2019), arXiv:1809.01116 [astro-ph.HE] .
  * Orsaria  _et al._ (2019) M. G. Orsaria, G. Malfatti, M. Mariani, I. F. Ranea-Sandoval, F. García, W. M. Spinella, G. A. Contrera, G. Lugones, and F. Weber, J. Phys. G 46, 073002 (2019), arXiv:1907.04654 [astro-ph.HE] .
  * Most  _et al._ (2020) E. R. Most, L. Jens Papenfort, V. Dexheimer, M. Hanauske, H. Stoecker, and L. Rezzolla, Eur. Phys. J. A 56, 59 (2020), arXiv:1910.13893 [astro-ph.HE] .
  * Yang  _et al._ (2021) W.-C. Yang, Y.-L. Ma, and Y.-L. Wu, Sci. China Phys. Mech. Astron. 64, 252011 (2021), arXiv:2011.03665 [nucl-th] .
  * Sekiguchi  _et al._ (2011) Y. Sekiguchi, K. Kiuchi, K. Kyutoku, and M. Shibata, Phys. Rev. Lett. 107, 211101 (2011), arXiv:1110.4442 [astro-ph.HE] .
  * Choi  _et al._ (2021) S. Choi, T. Miyatsu, M.-K. Cheoun, and K. Saito, Astrophys. J. 909, 156 (2021), arXiv:2011.14557 [nucl-th] .
  * Annala  _et al._ (2022) E. Annala, T. Gorda, E. Katerini, A. Kurkela, J. Nättilä, V. Paschalidis, and A. Vuorinen, Phys. Rev. X 12, 011058 (2022), arXiv:2105.05132 [astro-ph.HE] .
  * Raaijmakers  _et al._ (2021) G. Raaijmakers, S. K. Greif, K. Hebeler, T. Hinderer, S. Nissanke, A. Schwenk, T. E. Riley, A. L. Watts, J. M. Lattimer, and W. C. G. Ho, Astrophys. J. Lett. 918, L29 (2021), arXiv:2105.06981 [astro-ph.HE] .
  * Prakash  _et al._ (2021) A. Prakash, D. Radice, D. Logoteta, A. Perego, V. Nedora, I. Bombaci, R. Kashyap, S. Bernuzzi, and A. Endrizzi, Phys. Rev. D 104, 083029 (2021), arXiv:2106.07885 [astro-ph.HE] .
  * Kedia  _et al._ (2022) A. Kedia, H. I. Kim, I.-S. Suh, and G. J. Mathews, Phys. Rev. D 106, 103027 (2022), arXiv:2203.05461 [gr-qc] .
  * Fujimoto  _et al._ (2023) Y. Fujimoto, K. Fukushima, K. Hotokezaka, and K. Kyutoku, Phys. Rev. Lett. 130, 091404 (2023), arXiv:2205.03882 [astro-ph.HE] .
  * Guo  _et al._ (2023) L.-J. Guo, W.-C. Yang, Y.-L. Ma, and Y.-L. Wu, (2023), arXiv:2308.01770 [astro-ph.HE] .
  * Lenzi and Lugones (2012) C. H. Lenzi and G. Lugones, Astrophys. J. 759, 57 (2012), arXiv:1206.4108 [astro-ph.SR] .
  * Kojo (2021) T. Kojo, AAPPS Bull. 31, 11 (2021), arXiv:2011.10940 [nucl-th] .
  * Benic  _et al._ (2015) S. Benic, D. Blaschke, D. E. Alvarez-Castillo, T. Fischer, and S. Typel, Astron. Astrophys. 577, A40 (2015), arXiv:1411.2856 [astro-ph.HE] .
  * Minamikawa  _et al._ (2021a) T. Minamikawa, T. Kojo, and M. Harada, Phys. Rev. C 103, 045205 (2021a), arXiv:2011.13684 [nucl-th] .
  * Minamikawa  _et al._ (2021b) T. Minamikawa, T. Kojo, and M. Harada, Phys. Rev. C 104, 065201 (2021b), arXiv:2107.14545 [nucl-th] .
  * Gao  _et al._ (2022) B. Gao, T. Minamikawa, T. Kojo, and M. Harada, Phys. Rev. C 106, 065205 (2022), arXiv:2207.05970 [nucl-th] .
  * Kong  _et al._ (2023) Y. K. Kong, T. Minamikawa, and M. Harada, Phys. Rev. C 108, 055206 (2023), arXiv:2306.08140 [nucl-th] .
  * Gao  _et al._ (2024a) B. Gao, Y. Yan, and M. Harada, Phys. Rev. C 109, 065807 (2024a), arXiv:2404.04786 [nucl-th] .
  * DeTar and Kunihiro (1989) C. DeTar and T. Kunihiro, Phys. Rev. D 39, 2805 (1989).
  * Jido  _et al._ (2001) D. Jido, M. Oka, and A. Hosaka, Progress of Theoretical Physics 106, 873 (2001), https://academic.oup.com/ptp/article-pdf/106/5/873/5373808/106-5-873.pdf .
  * Aarts  _et al._ (2015) G. Aarts, C. Allton, S. Hands, B. Jäger, C. Praki, and J.-I. Skullerud, Phys. Rev. D 92, 014503 (2015), arXiv:1502.03603 [hep-lat] .
  * Aarts  _et al._ (2017) G. Aarts, C. Allton, D. De Boni, S. Hands, B. Jäger, C. Praki, and J.-I. Skullerud, JHEP 06, 034 (2017), arXiv:1703.09246 [hep-lat] .
  * Aarts  _et al._ (2019) G. Aarts, C. Allton, D. De Boni, and B. Jäger, Phys. Rev. D 99, 074503 (2019), arXiv:1812.07393 [hep-lat] .
  * Ma  _et al._ (2013) Y.-L. Ma, M. Harada, H. K. Lee, Y. Oh, B.-Y. Park, and M. Rho, Phys. Rev. D 88, 014016 (2013), [Erratum: Phys.Rev.D 88, 079904 (2013)], arXiv:1304.5638 [hep-ph] .
  * Ma  _et al._ (2014) Y.-L. Ma, M. Harada, H. K. Lee, Y. Oh, B.-Y. Park, and M. Rho, Phys. Rev. D 90, 034015 (2014), arXiv:1308.6476 [hep-ph] .
  * Nishihara and Harada (2015) H. Nishihara and M. Harada, Phys. Rev. D 92, 054022 (2015), arXiv:1506.07956 [hep-ph] .
  * Minamikawa  _et al._ (2023a) T. Minamikawa, B. Gao, T. Kojo, and M. Harada, (2023a), arXiv:2306.15564 [hep-ph] .
  * Hatsuda and Prakash (1989) T. Hatsuda and M. Prakash, Physics Letters B 224, 11 (1989).
  * Zschiesche  _et al._ (2007) D. Zschiesche, L. Tolos, J. Schaffner-Bielich, and R. D. Pisarski, Phys. Rev. C 75, 055202 (2007).
  * Dexheimer  _et al._ (2008) V. Dexheimer, S. Schramm, and D. Zschiesche, Phys. Rev. C 77, 025803 (2008).
  * Sasaki and Mishustin (2010) C. Sasaki and I. Mishustin, Phys. Rev. C 82, 035204 (2010).
  * Sasaki  _et al._ (2011) C. Sasaki, H. K. Lee, W.-G. Paeng, and M. Rho, Phys. Rev. D 84, 034011 (2011).
  * Gallas  _et al._ (2011) S. Gallas, F. Giacosa, and G. Pagliara, Nuclear Physics A 872, 13 (2011).
  * Steinheimer  _et al._ (2011) J. Steinheimer, S. Schramm, and H. Stöcker, Phys. Rev. C 84, 045208 (2011).
  * Paeng  _et al._ (2012) W.-G. Paeng, H. K. Lee, M. Rho, and C. Sasaki, Phys. Rev. D 85, 054022 (2012).
  * Dexheimer  _et al._ (2013) V. Dexheimer, J. Steinheimer, R. Negreiros, and S. Schramm, Phys. Rev. C 87, 015804 (2013).
  * Paeng  _et al._ (2013) W.-G. Paeng, H. K. Lee, M. Rho, and C. Sasaki, Phys. Rev. D 88, 105019 (2013).
  * Mukherjee  _et al._ (2017a) A. Mukherjee, J. Steinheimer, and S. Schramm, Phys. Rev. C 96, 025205 (2017a).
  * Suenaga (2018) D. Suenaga, Phys. Rev. C 97, 045203 (2018).
  * Takeda  _et al._ (2018) Y. Takeda, Y. Kim, and M. Harada, Phys. Rev. C 97, 065202 (2018).
  * Abuki  _et al._ (2018) H. Abuki, Y. Takeda, and M. Harada, EPJ Web Conf. 192, 00020 (2018).
  * Marczenko  _et al._ (2019a) M. Marczenko, D. Blaschke, K. Redlich, and C. Sasaki, Universe 5 (2019a), 10.3390/universe5080180.
  * Motohiro  _et al._ (2015) Y. Motohiro, Y. Kim, and M. Harada, Phys. Rev. C 92, 025201 (2015).
  * Yamazaki and Harada (2019) T. Yamazaki and M. Harada, Phys. Rev. C 100, 025205 (2019).
  * Mukherjee  _et al._ (2017b) A. Mukherjee, S. Schramm, J. Steinheimer, and V. Dexheimer, Astron. Astrophys. 608, A110 (2017b), arXiv:1706.09191 [nucl-th] .
  * Marczenko  _et al._ (2022) M. Marczenko, K. Redlich, and C. Sasaki, Phys. Rev. D 105, 103009 (2022), arXiv:2203.00269 [nucl-th] .
  * Minamikawa  _et al._ (2023b) T. Minamikawa, B. Gao, T. Kojo, and M. Harada, Symmetry 15, 745 (2023b), arXiv:2302.00825 [nucl-th] .
  * Marczenko  _et al._ (2019b) M. Marczenko, D. Blaschke, K. Redlich, and C. Sasaki, Universe 5, 180 (2019b), arXiv:1905.04974 [nucl-th] .
  * Gao  _et al._ (2024b) B. Gao, T. Kojo, and M. Harada, (2024b), arXiv:2403.18214 [hep-ph] .
  * Baym  _et al._ (2019) G. Baym, S. Furusawa, T. Hatsuda, T. Kojo, and H. Togashi, Astrophys. J. 885, 42 (2019), arXiv:1903.08963 [astro-ph.HE] .
  * Kobayashi and Maskawa (1970) M. Kobayashi and T. Maskawa, Prog. Theor. Phys. 44, 1422 (1970).
  * ’t Hooft (1986) G. ’t Hooft, Phys. Rept. 142, 357 (1986).
  * Hatsuda and Kunihiro (1994) T. Hatsuda and T. Kunihiro, Phys. Rept. 247, 221 (1994), arXiv:hep-ph/9401310 .
  * Zdunik and Haensel (2013) J. L. Zdunik and P. Haensel, Astron. Astrophys. 551, A61 (2013), arXiv:1211.1231 [astro-ph.SR] .
  * Alford  _et al._ (2013) M. G. Alford, S. Han, and M. Prakash, Phys. Rev. D 88, 083013 (2013), arXiv:1302.4732 [astro-ph.SR] .
  * Alford and Sedrakian (2017) M. G. Alford and A. Sedrakian, Phys. Rev. Lett. 119, 161104 (2017), arXiv:1706.01592 [astro-ph.HE] .
  * Glendenning and Kettner (2000) N. K. Glendenning and C. Kettner, Astron. Astrophys. 353, L9 (2000), arXiv:astro-ph/9807155 .