# Quark phases in neutron stars consistent with implications of NICER

Y. Yamamoto1<EMAIL_ADDRESS>N. Yasutake2 Th.A. Rijken31 1RIKEN Nishina
Center, 2-1 Hirosawa, Wako, Saitama 351-0198, Japan
2Department of Physics, Chiba Institute of Technology, 2-1-1 Shibazono
Narashino, Chiba 275-0023, Japan
3IMAPP, Radboud University, Nijmegen, The Netherlands

###### Abstract

The analyses for the NICER data imply $R_{2.0M_{\odot}}=12.41^{+1.00}_{-1.10}$
km and $R_{1.4M_{\odot}}=12.56^{+1.00}_{-1.07}$ km, indicating the lack of
significant variation of the radii from $1.4M_{\odot}$ to $2.0M_{\odot}$. This
feature cannot be reproduced by the hadronic matter due to the softening of
equation of state (EoS) by hyperon mixing, indicating the possible existence
of quark phases in neutron-star interiors. Two models are used for quark
phases: In the quark-hadron transition (QHT) model, quark deconfinement phase
transitions from a hadronic-matter EoS are taken into account so as to give
reasonable mass-radius ($MR$) curves by adjusting the quark-quark repulsions
and the density dependence of effective quark mass. In the quarkyonic model,
the degrees of freedom inside the Fermi sea are treated as quarks and neutrons
exist at the surface of the Fermi sea, where $MR$ curves are controlled mainly
by the thickness of neutron Fermi layer. The QHT and quarkyonic EoSs can be
adjusted so as to reproduce radii, tidal deformabilities, pressure and central
densities inferred from the NICER analysis better than the nucleonic matter
EoS, demonstrating the clear impacts of quark phases. Then, the maximum mass
for the quakyonic-matter EoS is considerably larger than that for the QHT-
matter EoS.

###### pacs:

21.30.Cb, 21.30.Fe, 21.65.+f, 21.80.+a, 12.39.Jh, 25.75.Nq, 26.60.+c

##  I Introduction

In studies of neutron stars (NS), the fundamental role is played by the
equation of state (EoS) for neutron star matter. The massive neutron stars
with masses over $2M_{\odot}$ have been reliably established by the
observations of NSs J1614$-$2230 Demorest10 , J0348+0432 Antoniadis13 ,
J0740+6620 Cromartie2020  and J0952-0607 Romani2022 . The radius information
of NSs have been obtained for the massive NS PSR J0740+6620 with $2M_{\odot}$
and $1.4M_{\odot}$ NSs, shown as $R_{2M_{\odot}}$ and $R_{1.4M_{\odot}}$, from
the analyses for the X-ray data taken by the Neutron Star Interior Composition
Explorer (NICER) and the X-ray Multi-Mirror (XMM-Newton) observatory. The
analysis of Miller et al. gives $R_{2.08M_{\odot}}=12.35\pm 0.75$ km and
$R_{1.4M_{\odot}}=12.45\pm 0.65$ km Miller2021 . The analysis of Riley gives
$R_{2.08M_{\odot}}=12.39^{+1.30}_{-0.98}$ km and
$R_{1.4M_{\odot}}=12.33^{+0.76}_{-0.81}$ km Riley2021 . Legred et al.
investigate these measurement’s implications for the EoSs, employing a
nonparametric EoS model based on Gaussian processes and combining information
from other X-ray and gravitational wave observations Legred2021 .

The purpose of this paper is to demonstrate that the radius information of
massive NSs give the important constraints for the neutron-star EoSs. In our
EoS analysis, the following neutron-star radii are adopted as critical values
to be reproduced:

$\displaystyle R_{2.0M_{\odot}}=12.41^{+1.00}_{-1.10}\ \ {\rm km}$
$\displaystyle R_{1.4M_{\odot}}=12.56^{+1.00}_{-1.07}\ \ {\rm km}$ (1)

with maximum mass $M_{max}/M_{\odot}=2.21^{+0.31}_{-0.21}$, being given by the
analysis by Legred et al.Legred2021 . The median values of $R_{2M_{\odot}}$
and $R_{1.4M_{\odot}}$ in the above three references Miller2021 Riley2021
Legred2021  are only a few hundred meters apart from each other. We set the
fitting accuracy to a few hundred meters in our analysis for $R_{2M_{\odot}}$
and $R_{1.4M_{\odot}}$. Then, the EoS obtained from our analysis are not
changed, even if the set of $R_{2M_{\odot}}$ and $R_{1.4M_{\odot}}$ in
Miller2021  or Riley2021  is used as the criterion instead of Eq.(1) or all
three sets in Miller2021 Riley2021 Legred2021  are used. The key feature found
commonly in the three sets is the small variation of radii from $1.4M_{\odot}$
to $2M_{\odot}$, namely $R_{2M_{\odot}}\approx R_{1.4M_{\odot}}$. The reason
why the result in Legred2021  is used in our present analysis is because they
present the inferred values of maximum masses, radii, tidal deformabilities,
pressure and central densities obtained from their analysis. These quantities
can be compared with our corresponding results, by which the features of our
EoSs are revealed in detail.

The hyperon mixing in neutron-star matter brings about a remarkable softening
of the EoS and a maximum mass is reduced to a value far less than
$2M_{\odot}$. The EoS softening is caused by changing of high-momentum
neutrons at Fermi surfaces to low-momentum hyperons via strangeness non-
conserving weak interactions overcoming rest masses of hyperons. In order to
derive EoSs for massive NSs, it is necessary to solve this “hyperon puzzle in
neutron stars”. There have been proposed possible mechanisms: (i) more
repulsive hyperon-hyperon interactions in relativistic mean field (RMF) models
driven by vector mesons exchanges Bednarek ; Weiss ; Oertel ; Maslov , (ii)
repulsive hyperonic three-body forces NYT ; Vidana11 ; YFYR14 ; YFYR16 ;
YTTFYR17 ; Lonardoni ; Logoteta ; Cerstung , (iii) appearance of other
hadronic degrees of freedom, such as $\Delta$ isobars Drago  or meson
condensates Kaplan ; Brown ; Thorsson ; Lee ; GleSch , (iv) existence of quark
phases in high-density regions Baldo2006 ; Ozel ; Weissenborn ; Klahn ;
Bonanno ; Lastowiecki2012 ; Shahrbaf1 ; Shahrbaf2 ; Otto2020 ; KBH2022 ;
YYR2022 . It should be noted that the criterion for NS radii Eq.(1) is
stricter than the condition of $M_{max}>2M_{\odot}$ only to solve the ”puzzle”
and the above mechanisms are needed to be re-investigated under this stricter
condition.

One of the approaches belonging to (ii) is to assume that three-nucleon
repulsions (TNR) APR98  work universally among every kind of baryons as three-
baryon repulsions (TBR) NYT . In YFYR14 ; YFYR16 ; YTTFYR17 , the multi-
pomeron exchange potential (MPP) was introduced as a model of universal
repulsions among three and four baryons on the basis of the extended soft core
(ESC) baryon-baryon interaction model developed by two of the authors (T.R.
and Y.Y.) and M.M. Nagels ESC16I ; ESC16II ; ESC16III . In the case of this
special modeling for hyperonic three-body repulsions, the EoS softening by
hyperon mixing is not completely recovered by the above universal repulsions,
and the maximum masses become not so large even if universal many-body
repulsions increase. As a result, the maximum masses for hyperonic-matter EoS
cannot be over $2M_{\odot}$, as found in YFYR14 ; YFYR16 ; YTTFYR17 : It is
difficult that criterion Eq.(1) is realized by this modeling of hadronic-
matter EoSs. A simple way to avoid the strong softening of EoS by hyperon
mixing is to assume $\Lambda NN$ repulsions stronger than $NNN$ repulsions
with neglect of $\Sigma^{-}$ mixing Lonardoni .

In this paper, we focus on the mechanism (iv). It is possible to solve the
“hyperon puzzle” by taking account of quark deconfinement phase transitions
from a hadronic-matter EoS to a sufficiently stiff quark-matter EoS in the
neutron-star interiors, namely by studying hybrid stars having quark matter in
their cores, where repulsive effects in quark phases are needed to result in
massive stars over $2M_{\odot}$. In the Nambu-Jona-Lasinio (NJL) model, for
instance, repulsions to stiffen EoSs are given by vector interactions. Then,
it is known well that quark-hadron phase transitions should be crossover or at
most weak first-order, because strong first-order transitions soften EoSs
remarkably in order to obtain stiff EoSs. In KBH2022 , they derived the new
EoS within the quark-hadron crossover (QHC) framework (3-windows model) so as
to reproduce $R_{2.1M_{\odot}}\approx R_{1.4M_{\odot}}\approx 12.4{\rm km}$.
Here, the small variation of radii indicates that the pressure grows rapidly
while changes in energy density are modest, producing a peak in the speed of
sound KBH2022 . In their QHC framework, the EoSs in the quark-hadron mixed
region of $1.5\rho_{0}\sim 3.5\rho_{0}$, playing a decisive role for the
resulting $MR$ curves, are given by the interpolating functions
phenomenologically. Then, it is meaningful to study the other modeling for
phase transitions in which the mixed regions are modeled explicitly. We
investigate how this criterion Eq.(1) can be realized in the case of using the
EoS derived from our quark-hadron transition (QHT) model for neutron-star
matter in the Bruecner-Hartree-Fock (BHF) framework YYR2022 , being different
from their 3-windows model. Here, the quark-matter EoS is derived from the
two-body quark-quark ($QQ$) potentials, in which all parameters are on the
physical backgrounds with no room for arbitrarily changing: They are composed
of meson-exchange quark-quark potentials derived by unfolding of the baryon-
baryon meson-exchanges, and instanton-exchange, one-gluon-exchange and multi-
pomeron exchange potentials. Then, baryonic matter and quark matter are
treated in the common BHF framework, where quark-hadron transitions are
treated on the basis of the Maxwell condition. In this paper, it is shown that
the criterion Eq.(1) can be realized by our QHT model for neutron-star matter,
as well as the QHC model KBH2022 , by adjusting the $QQ$ repulsion to be
strong enough and the quark-hadron transition density to be about $2\rho_{0}$.

In our QHT model the BHF framework is used for deriving the quark-matter EoS,
which is not popular. Our treatments for quark-hadron phase transitions is the
same as that in Shahrbaf2  where the NJL model is adopted for quark matter
under the mean field approximation. In spite of the difference between quark-
matter models, their obtained $MR$ curves are similar to ours in YYR2022 .
Therefore, it is considered that the same conclusions can be derived also by
using their QHT model instead of ours.

Another type of quark phase in neutron-star interiors is given by the
quarkyonic matter MP2007 ; HMP2008 ; MR2019 ; HMLCP2019 ; DHJ2020 ; ZL2020 ;
MHPC2021 ; Cao2022 , where the degrees of freedom inside the Fermi sea are
treated as quarks and nucleons exist at the surface of the Fermi sea. The
transition from hadronic-matter phase to the quarkyonic-matter phase is
considered to be in second-order. In the quarkyonic matter, the existence of
free quarks inside the Fermi sea gives nucleons extra kinetic energy by
pushing them to higher momenta, leading to increasing pressure. This mechanism
to realize the criterion Eq.(1) is completely different from the QHT matter in
which the essential roles for EoS stiffening are played by the $QQ$
repulsions. Then, it is valuable to study the characteristic differences
between neutron-star mass-radius ($MR$) curves obtained from the QHT-matter
EoS and quarkyonic-matter EoS.

This paper is organized as follows: In Sect.II, the hadronic-matter EoS
(II-A), the quark-matter EoS (II-B) and the quarkyonic-matter EoS (II-C) are
formulated on the basis of our previous works, where the BHF frameworks with
our $QQ$ potentials are adopted both for baryonic matter and quark
(quarkyonic) matter. Transitions from hadron phases to quark matter
(quakyonic) phases are explained. In Sect.III-A, the calculated results are
shown for pressures, energy densities and sound velocities. In III-B, the $MR$
curves of hybrid stars are obtained by solving the Tolmann-Oppenheimer-Volkoff
(TOV) equation. In III-C, the obtained values of maximum masses, radii, tidal
deformabilities, pressure and central densities are compared with those
inferred from the NICER-data analysis. The conclusion of this paper is given
in Sect.IV.

##  II Models of neutron-star matter

###  II.1 hadronic matter

The hadronic matter is defined as $\beta$-stable hyperonic nuclear matter
including leptons, composed of $n$, $p^{+}$, $\Lambda$, $\Sigma^{-}$, $e^{-}$,
$\mu^{-}$. We recapitulate here the hadronic-matter EoS. In the BHF framework,
the EoS is derived with use of the ESC baryon-baryon ($B\\!B$) interaction
model YFYR14 ; YFYR16 ; YTTFYR17 .

As is well known, the nuclear-matter EoS is stiff enough to assure neutron-
star masses over $2M_{\odot}$, if the strong three-nucleon repulsion (TNR) is
taken into account. However, there appears a remarkable softening of EoS by
inclusion of exotic degrees of freedom such as hyperon mixing. One of the
ideas to avoid this “hyperon puzzle” is to assume that the many-body
repulsions work universally for every kind of baryons NYT . In YFYR14 ; YFYR16
; YTTFYR17 , the multi-pomeron exchange potential MPP was introduced as a
model of universal repulsions among three and four baryons. This was inspired
by the multi-reggeon model to describe CERN-ISR pp-data KM74 . The ESC work is
mentioned in ESC16I ; ESC16II ; ESC16III .

In YTTFYR17  they proposed three versions of MPP (MPa, MPa+, MPb), where MPa
and MPa+ (MPb) include the three- and four-body (only three-body) repulsions.
Their strengths are determined by analyzing the nucleus-nucleus scattering
using the G-matrix folding model under the conditions that the saturation
parameters are reproduced reasonably. The EoSs for MPa and MPa+ are stiffer
than that for MPb, and maximum masses and radii of neutron stars obtained from
MPa, MPa+ are larger than those from MPb. The important criterion for
repulsive parts is the resulting neutron-star radii $R$ for masses of
$1.4M_{\odot}$: In the case of using MPb, we obtain $R_{1.4M_{\odot}}\approx
12.4$ km similar to the value in the criterion Eq.(1). On the other hand, we
have $R_{1.4M_{\odot}}\approx 13.3$ (13.6) km in the case of MPa (MPa+). In
this paper, we adopt MPb as three-baryon repulsion: Our nuclear interactions
are composed of two-body part $V_{BB}$ and three-body part $V_{BBB}$, where
$V_{BB}$ and $V_{BBB}$ are given by ESC and MPb, respectively. It is
worthwhile to say that the three-nucleon repulsion in MPb is stronger than the
corresponding one (UIX) in the standard model by APR APR98  giving rise to
$R_{1.4M_{\odot}}\approx 11.6$ km Togashi1 .

$BB$ G-matrix interactions ${\cal G}_{BB}$ are derived from $BB$ bare
interactions $V_{BB}$ or $V_{BB}+V_{BBB}$ YFYR14 . They are given for each
$(BB^{\prime},T,S,P)$ state, $T$, $S$ and $P$ being isospin, spin and parity
in a two-body state, respectively, and represented as ${\cal
G}_{BB^{\prime}}^{TSP}$. The G-matrix interactions derived from $V_{BB}$ and
$V_{BB}+V_{BBB}$ are called B1 and B2, respectively. In the quarkyonic model,
we need only the neutron-neutron sectors, ${\cal G}_{nn}^{SP}$.

A single baryon potential is given by

$\displaystyle U_{B}(k)$ $\displaystyle=$
$\displaystyle\sum_{B^{\prime}=n,p,\Lambda,\Sigma^{-}}U_{B}^{(B^{\prime})}(k)$
$\displaystyle=$
$\displaystyle\sum_{B^{\prime}=n,p,\Lambda,\Sigma^{-}}\sum_{k^{\prime}<k_{F}^{(B^{\prime})}}\langle
kk^{\prime}|{\cal G}_{BB^{\prime}}|kk^{\prime}\rangle$

with $B=n,p,\Lambda,\Sigma^{-}$. Here, $\langle kk^{\prime}|{\cal
G}_{BB^{\prime}}|kk^{\prime}\rangle$ is a $BB^{\prime}$ G-matrix element in
momentum space, being derived from $V_{BB}$ or ($V_{BB}$+$V_{BBB}$), and
$k_{F}^{(B)}$ is the Fermi momentum of baryon $B$. In this expression, spin
and isospin quantum numbers are implicit.

The baryon energy density is given by

$\displaystyle\varepsilon_{B}$ $\displaystyle=$
$\displaystyle\tau_{B}+\upsilon_{B}$ $\displaystyle=$ $\displaystyle
g_{s}\int_{0}^{k_{F}^{(B)}}\frac{d^{3}k}{(2\pi)^{3}}\left\\{\sqrt{\hbar^{2}k^{2}+M_{B}^{2}}+\frac{1}{2}U_{B}(k)\right\\}\
,$

where $\tau_{B}$ and $\upsilon_{B}$ are kinetic and potential parts of the
energy density.

In $\beta$-stable hadronic matter composed of $n$, $p$, $e^{-}$, $\mu^{-}$,
$\Lambda$ and $\Sigma^{-}$, equilibrium conditions are given as

(1) chemical equilibrium conditions,

$\displaystyle\mu_{n}=\mu_{p}+\mu_{e}$ (4) $\displaystyle\mu_{\mu}=\mu_{e}$
(5) $\displaystyle\mu_{\Lambda}=\mu_{n}$ (6)
$\displaystyle\mu_{\Sigma^{-}}=\mu_{n}+\mu_{e}$ (7)

(2) charge neutrality,

$\displaystyle\rho_{p}=\rho_{e}+\rho_{\mu}+\rho_{\Sigma^{-}}$ (8)

(3) baryon number conservation,

$\displaystyle\rho=\rho_{n}+\rho_{p}+\rho_{\Lambda}+\rho_{\Sigma^{-}}\ .$ (9)

Expressions for $\beta$-stable nucleonic matter composed of $n$, $p$, $e^{-}$
and $\mu^{-}$ are obtained by omitting hyperon sectors from the above
expressions for $\beta$-stable baryonic matter.

###  II.2 Quark-Hadron transition model

In our treatment of quark matter, the BHF framework is adopted on the basis of
two-body $QQ$ potentials YYR2022 . Here, correlations induced by bare $QQ$
potentials are renormalized into coordinate-space G-matrix interactions, being
considered as effective $QQ$ interactions used in quark-matter calculations.

Our bare $QQ$ interaction is given by

$\displaystyle V_{QQ}$ $\displaystyle=$ $\displaystyle
V_{EME}+V_{INS}+V_{OGE}+V_{MPP}$ (10)

where $V_{EME}$, $V_{INS}$, $V_{OGE}$ and $V_{MPP}$ are the extended meson-
exchange potential, the instanton-exchange potential, the one-gluon exchange
potential and the multi-pomeron exchange potential, respectively. Parameters
in our $QQ$ potential are chosen so as to be consistent with physical
observables. The $V_{EME}$ $QQ$ potential is derived from the ESC $BB$
potential so that the $QQM$ couplings are related to the $BBM$ couplings
through folding procedures with Gaussian baryonic quark wave functions. In the
construction of the relation between $BBM$ and $QQM$ couplings, the
requirement that the coefficients of the $1/M^{2}$ expansion should match is
based on Lorentz invariance, which fixes the QQM couplings and also determines
the (few) extra vertices at the quark level ESC16I . Then, the $V_{EME}$ $QQ$
potential is basically of the same functional expression as the ESC $BB$
potential. Strongly repulsive components in ESC $BB$ potentials are described
mainly by vector-meson and pomeron exchanges between baryons. This feature
persists in the $V_{EME}$ $QQ$ potential, which includes the strongly
repulsive components originated from vector-meson and pomeron exchanges
between quarks. Similarly the multi-pomeron exchange potentials among quarks,
$V_{MPP}$, are derived from the corresponding ones among baryons, giving
repulsive contributions. Contributions from $V_{INS}$ and $V_{OGE}$ in average
are attractive and repulsive, respectively. The strength of $V_{OGE}$ is
determined by the quark-gluon coupling constant $\alpha_{S}$. In YYR2022
$\alpha_{S}$ is chosen as 0.25, that is $V_{OGE}(\alpha_{S}=0.25)$, and the
three sets are defined as follows: Q0 : $V_{EME}$, Q1 :
$V_{EME}+V_{INS}+V_{OGE}(\alpha_{S}=0.25)$ Q2 :
$V_{EME}+V_{MPP}+V_{INS}+V_{OGE}(\alpha_{S}=0.25)$.

In our QHT model for neutron-star matter, quark-hadron phase transitions occur
at crossing points of hadron pressure $P_{H}(\mu)$ and quark pressure
$P_{Q}(\mu)$ being a function of chemical potential $\mu$. Positions of
crossing points, giving quark-hadron transition densities, are controlled by
parameters $\rho_{c}$ and $\gamma$ included in our density-dependent quark
mass

$\displaystyle
M_{Q}^{*}(\rho_{Q})=\frac{M_{0}}{1+\exp[\gamma(\rho_{Q}-\rho_{c})]}+m_{0}+C$
(11)

with $C=M_{0}-M_{0}/[1+\exp(-\gamma\rho_{c})]$ assuring
$M_{Q}^{*}(0)=M_{0}+m_{0}$, where $\rho_{Q}$ is number density of quark
matter, and $M_{0}$ and $m_{0}$ are taken as 300 (360) MeV and 5 (140) MeV for
$u$ and $d$ ($s$) quarks. Here, the effective quark mass $M_{Q}^{*}(\rho_{Q})$
should be used together with
$B(\rho_{Q})=M_{Q}^{*}(0)-M_{Q}^{*}(\rho_{Q})+B_{0}$, meaning the energy-
density difference between the perturbative vacuum and the true vacuum. A
constant term $B_{0}$ is added for fine tuning of an onset density. In YYR2022
, the values of ($\rho_{c}$, $\gamma$) without $B_{0}$ are given for each set
of Q0, Q1 and Q2.

Let us focus on the typical result for Q2+H1 in YYR2022 . The $QQ$ interaction
Q2 is the most repulsive among Q0, Q1 and Q2. The $BB$ interaction H1 consists
of ESC and MPb, and results in the reasonable value of $R_{1.4M_{\odot}}$. In
this case of Q2+H1, we obtain the maximum mass of $2.25M_{\odot}$ and the
reasonable value of $R_{1.4M_{\odot}}=12.5$ km, in which the quark-hadron
transition occurs at density of 3.5$\rho_{0}$. Then, we have
$R_{2.0M_{\odot}}=12.0$ km, being rather smaller than 12.4 km in the criterion
Eq.(1). In order to reproduce a larger value of $R_{2.0M_{\odot}}\approx 12.4$
km, we make $V_{OGE}$ more repulsive by taking larger values of $\alpha_{S}=$
0.36 and 0.49. It is not suitable for such a purpose to strengthen the
$V_{MPP}$ repulsion, because $V_{MPP}$ is essentially of three-body
interaction and the contributions in low-density region are small. On the
other hand, $V_{OGE}$ is of two-body interaction, and its repulsive
contributions are not small even in low density region, being important for a
large value of $R_{2.0M_{\odot}}$. Another condition to make
$R_{2.0M_{\odot}}$ larger is to lower quark-hadron transition densities by
adjusting the parameters ($\rho_{c}$,$\gamma$,$B_{0}$) included in the
density-dependent quark mass Eq.(11).

We define newly the following three sets with the fixed value of $\gamma$=1.2

Q2 : $V_{EME}+V_{MPP}+V_{INS}+V_{OGE}(\alpha_{S}=0.25)$

with $\rho_{c}=6.9\rho_{0}$ and $B_{0}=$8.5

Q3 : $V_{EME}+V_{MPP}+V_{INS}+V_{OGE}(\alpha_{S}=0.36)$

with $\rho_{c}=6.9\rho_{0}$ and $B_{0}=$7.5

Q4 : $V_{EME}+V_{MPP}+V_{INS}+V_{OGE}(\alpha_{S}=0.69)$

with $\rho_{c}=7.5\rho_{0}$ and $B_{0}=$10.0

where the values of $\rho_{c}$ and $B_{0}$ for each set are chosen so as to
give quark-hadron transition densities of $\sim 2\rho_{0}$.

G-matrix interactions ${\cal G}_{qq^{\prime}}$ with $q,q^{\prime}=u,d,s$ are
derived from the above bare $QQ$ interactions. They are given for each
$(qq^{\prime},T,S,P)$ state, $T$, $S$ and $P$ being isospin, spin and parity
in a two-body state, respectively, and represented as ${\cal
G}_{qq^{\prime}}^{TSP}$. Hereafter, Q2, Q3 and Q4 mean the naming of
corresponding $QQ$ G-matrix interactions, not only of bare $QQ$ interactions.
The $QQ$ G-matrix interactions are used also in the quarkyonic matter
calculations.

A single quark potential is given by

$\displaystyle U_{q}(k)$ $\displaystyle=$
$\displaystyle\sum_{q^{\prime}=u,d,s}U_{q}^{(q^{\prime})}(k)=\sum_{q^{\prime}=u,d,s}\sum_{k^{\prime}<k_{F}^{q^{\prime}}}\langle
kk^{\prime}|{\cal G}_{qq^{\prime}}|kk^{\prime}\rangle$

with $q=u,d,s$, where $k_{F}^{q}$ is the Fermi momentum of quark $q$. Spin and
isospin quantum numbers are implicit.

The quark energy density is given by

$\displaystyle\varepsilon_{q}$ $\displaystyle=$ $\displaystyle
g_{s}N_{c}\sum_{q=u,d,s}\int_{0}^{k_{Fq}}\frac{d^{3}k}{(2\pi)^{3}}$
$\displaystyle\left\\{\sqrt{\hbar^{2}k^{2}+M_{q}^{2}}+\frac{1}{2}U_{q}(k)\right\\}\
.$

Fermion spin and quark color degeneracies give rise to $g_{s}=2$ and
$N_{c}=3$.

Figure 1: (Color online) Potential energies per particle $U/A$ as a function
of the baryon number density $\rho_{B}$ in the case of
$\rho_{u}=\rho_{d}=\rho_{s}$. The short-dashed, long-dashed and solid curves
are obtained by using Q2, Q3 and Q4, respectively.

In order to demonstrate the features of our $QQ$ interactions (Q2,Q3,Q4), we
show the potential energy per particle $U/A$ as a function of the baryon
number density $\rho_{B}=\frac{1}{3}\rho_{Q}$ in the case of taking
$\rho_{u}=\rho_{d}=\rho_{s}$. In Fig.1, the short-dashed, long-dashed and
solid curves are obtained by using Q2, Q3 and Q4, respectively. The repulsions
are found to be strong in the order of Q4, Q3, Q2. This difference of
repulsions among Q4, Q3 and Q2 comes from the different values of $\alpha_{S}$
included in $V_{OGE}$. In the figure, it should be noted that the difference
is not small even in the low-density region.

In the EoS of $\beta$-stable quark matter composed of $u$, $d$, $s$, $e^{-}$,
the equilibrium conditions are given as

(1) chemical equilibrium conditions,

$\displaystyle\mu_{d}=\mu_{s}=\mu_{u}+\mu_{e}$ (14)

(2) charge neutrality,

$\displaystyle 0=\frac{1}{3}(2\rho_{u}-\rho_{d}-\rho_{s})-\rho_{e}$ (15)

(3) baryon number conservation,

$\displaystyle\rho_{B}=\frac{1}{3}(\rho_{u}+\rho_{d}+\rho_{s})=\frac{1}{3}\rho_{Q}\
.$ (16)

In order to construct the hybrid EoS including a transition from hadronic
phase to quark phase, we use the replacement interpolation method Shahrbaf2
YYR2022 , being a simple modification of the Maxwell and the Glendenning
(Gibbs) constructions Glendenning . The EoSs of hadronic and quark phases and
that of mixed phase are described with the relations between pressures and
chemical potentials $P_{H}(\mu)$, $P_{Q}(\mu)$ and $P_{M}(\mu)$, respectively.
The critical chemical potential $\mu_{c}$ for the transition from the hadronic
phase to the quark phase is obtained from the Maxwell condition

$\displaystyle P_{Q}(\mu_{c})=P_{H}(\mu_{c})=P_{c}\ .$ (17)

The pressure of the mixed phase is represented by a polynomial ansatz. The
matching densities $\rho_{H}$ and $\rho_{Q}$ are obtained with use of
$\rho(\mu)=dP(\mu)/d\mu$.

###  II.3 quarkyonic matter

In the BHF framework, we derive the EoS of quarkyonic matter composed of
neutrons and quarks with flavor $q=u,d$ in the simplest form by McLerran and
Reddy MR2019 . In the chargeless 2-flavor quarkyonic matter, strongly
interacting quarks near the Fermi sea form interacting neutrons, and the
remaining d and u quarks fill the lowest momenta up to $k_{Fu}$ and $k_{Fd}$,
respectively. The quark mass is taken to be $M_{q}=M_{n}/3$ constantly,
$M_{n}$ being the neutron mass. In calculations of quarkyonic matter, we use
B1 ($V_{nn}$) and B2 ($V_{nn}$+$V_{nnn}$) for nuclear interactions, and Q0 for
$QQ$ interactions for simplicity.

The total baryon number density is given by

$\displaystyle\rho_{B}$ $\displaystyle=$
$\displaystyle\rho_{n}+\frac{N_{c}}{3}(\rho_{u}+\rho_{d})$ (18)
$\displaystyle=$
$\displaystyle\frac{g_{s}}{6\pi^{2}}\left[k_{Fn}^{3}-k_{0n}^{3}+\frac{N_{c}}{3}(k_{Fu}^{3}+k_{Fd}^{3})\right]\
,$

where $k_{Fn}$, $k_{Fu}$ and $k_{Fd}$ are the Fermi momenta of neutrons and u
and d quarks, respectively. Fermion spin and quark color degeneracies give
rise to $g_{s}=2$ and $N_{c}=3$. Neutrons are restricted near the Fermi
surface by $k_{0n}$, being assumed as

$\displaystyle k_{0n}=k_{Fn}-\Delta_{qyc}$
$\displaystyle\Delta_{qyc}=\frac{\Lambda^{3}}{\hbar
c^{3}k_{Fn}^{2}}+\kappa\frac{\Lambda}{N_{c}^{2}\hbar c}\ ,$ (19)

where $\Delta_{qyc}$ for the thickness of Fermi layer includes the two
parameters $\Lambda$ and $\kappa$. In this work, we take the fixed value of
$\kappa=0.3$.

Then, $k_{Fd}$ and $k_{Fu}$ are related to $k_{0n}$ by
$k_{Fd}=\frac{1}{N_{c}}k_{0n}$ and $k_{Fu}=2^{-1/3}k_{Fd}$.

A single neutron potential is given by

$\displaystyle U_{n}(k)$ $\displaystyle=$
$\displaystyle\sum_{k_{0n}<k^{\prime}<k_{Fn}}\langle kk^{\prime}|{\cal
G}_{nn}|kk^{\prime}\rangle$ (20)

with $nn$ G-matrix interactions ${\cal G}_{nn}$.

The neutron energy density is given by

$\displaystyle\varepsilon_{n}$ $\displaystyle=$
$\displaystyle\tau_{n}+\upsilon_{n}$ $\displaystyle=$ $\displaystyle
g_{s}\int_{k_{0n}}^{k_{Fn}}\frac{d^{3}k}{(2\pi)^{3}}\left\\{\sqrt{\hbar^{2}k^{2}+M_{n}^{2}}+\frac{1}{2}U_{n}(k)\right\\}\
.$

Additionally, another form of the neutron potential energy density is defined
as

$\displaystyle{\bar{\upsilon}_{n}}=g_{s}\int_{0}^{k_{n}}\frac{d^{3}k}{(2\pi)^{3}}\left\\{\frac{1}{2}U_{n}(k)\right\\}\
,$ (22)

which is used in MR2019  instead of $\upsilon_{n}$.

Single quark potentials for $q=u,d$ are given by

$\displaystyle U_{q}(k)$ $\displaystyle=$
$\displaystyle\sum_{q^{\prime}=u,d}U_{q}^{(q^{\prime})}(k)$ (23)
$\displaystyle=\sum_{q^{\prime}=u,d}\sum_{k^{\prime}<k_{Fq}}\langle
kk^{\prime}|{\cal G}_{qq^{\prime}}|kk^{\prime}\rangle$ $\displaystyle
U_{q}^{(n)}(k)$ $\displaystyle=$
$\displaystyle\sum_{k_{0n}<k^{\prime}<k_{Fn}}\langle kk^{\prime}|{\cal
G}_{qn}|kk^{\prime}\rangle$ (24)

with G-matrix interactions ${\cal G}_{qq^{\prime}}$ and ${\cal G}_{qn}$. Here,
${\cal G}_{qn}$ is the quark-neutron ($Qn$) interactions: We assume the simple
model in which the potentials ${\cal G}_{qq^{\prime}}$ are folded into the
potentials ${\cal G}_{qn}$ with Gaussian baryonic quark wave functions. In
Eqs.(20)(23)(24) spin quantum numbers are implicit.

The quark energy density is given by

$\displaystyle\varepsilon_{q}$ $\displaystyle=$ $\displaystyle
g_{s}N_{c}\sum_{q=u,d}\int_{0}^{k_{Fq}}\frac{d^{3}k}{(2\pi)^{3}}$
$\displaystyle\left\\{\sqrt{\hbar^{2}k^{2}+M_{q}^{2}}+\frac{1}{2}U_{q}(k)+U_{qn}(k)\right\\}\
,$

where values of $k_{Fq}$ are determined by

$\displaystyle N_{c}k_{Fq}=k_{0n}\ .$ (26)

Thus, our total energy density is given by

$\displaystyle\varepsilon=\varepsilon_{n}+\varepsilon_{d}+\varepsilon_{u}\ .$
(27)

The chemical potential $\mu_{i}$ ($i=n,d,u$) and pressure $P$ are expressed as

$\displaystyle\mu_{i}=\frac{\partial\varepsilon_{i}}{\partial n_{i}}\ ,$ (28)
$\displaystyle P=\sum_{i=n,d,u}\mu_{i}n_{i}-\varepsilon\ ,$ (29)

where $\frac{\partial\varepsilon_{i}}{\partial
n_{i}}=\frac{\partial\varepsilon_{i}}{\partial n_{B}}\frac{\partial
n_{B}}{\partial n_{i}}$ .

In our model, the phase transition from $\beta$-stable nucleonic matter to the
quarkyonic matter occurs in second-order, resulting in the hybrid EoS
including hadronic and quarkyonic EoSs. Then, the transition densities are
controlled mainly by the parameter $\Lambda$: In this work, we choose the
three values of $\Lambda$=380, 350 and 320 MeV with the fixed value of
$\kappa=0.3$. The transition densities for these values are $0.28\sim 0.38$
fm-3 ($0.28\sim 0.36$ fm-3) in the case of using B1 (B2) for nuclear
interactions. Hereafter, when a value of $\Lambda$=380 MeV is used, for
instance, it is denoted as $\Lambda 380$.

##  III Results and discussion

###  III.1 EoS

Figure 2: (Color online) Pressures $P$ as a function of baryonic number
density $\rho_{B}$. The dot-dashed (dotted) curve is for $\beta$-stable
nucleonic (hadronic) matter. Upper thin (thick) solid curves are pressures in
the quarkyonic matter for $\Lambda 350$ and $\Lambda 320$ ($\Lambda 380$) with
B1. Lower thin (thick) short-dashed curves are for the QHT matter with Q2 and
Q3 (Q4). Figure 3: (Color online) Pressures $P$ as a function of the energy
density $\varepsilon$. The dot-dashed (dotted) curves are for $\beta$-stable
nucleonic (hadronic) matter. Thin (thick) solid curves show pressures in
quarkyonic phases for $\Lambda 350$ and $\Lambda 320$ ($\Lambda 380$) with B1.
The short-dashed curve is for the QHT model with Q4. Figure 4: (Color online)
The square of the sound speed $c_{s}^{2}$ in units of $c^{2}$ as a function of
baryonic number density $\rho_{B}$. The dot-dashed (dotted) curve is that in
$\beta$-stable nucleonic (hadronic) matter. Solid curves are pressures in
quarkyonic matter for $\Lambda 380$, $\Lambda 350$ and $\Lambda 320$ with B1.
The dashed curve is for the QHT matter with Q4. Figure 5: (Color online) In
the left panel, solid curves are pressures $P$ in quarkyonic phases for as a
function of baryonic number density $\rho_{B}$ for $\Lambda 380$ in the cases
of using B1 and B2, and short-dashed (dashed) curves are partial pressures of
neutrons (quarks) in respective cases. The dot-dashed curve is for
$\beta$-stable nucleonic matter. In the right panel, solid (dot-dashed) curves
are neutron chemical potentials $\mu_{n}$ in quarkyonic ($\beta$-stable
nucleonic) phases as a function of $\rho_{B}$ for $\Lambda 380$ in the cases
of using B1 and B2. The dot-dashed curve gives neutron chemical potential in
$\beta$-stable nucleonic matter.

In Fig.2, pressures $P$ are drawn as a function of baryonic number density
$\rho_{B}$. The dot-dashed curve is for the $\beta$-stable nucleonic-matter
EoS, and the dotted one is for the $\beta$-stable hadronic-matter EoS with
hyperon mixing. The latter is substantially below the former, demonstrating
the EoS softening by hyperon mixing. Thin (thick) solid curves in the upper
side are pressures in the quarkyonic matter for $\Lambda 350$ and $\Lambda
320$ ($\Lambda 380$) with use of B1 for nuclear interactions. At the crossing
points with the dot-dashed curve in the low-density side, there occur second-
order transitions from $\beta$-stable nucleonic to quarkyonic phases: The
transition densities $\rho_{t}$ are 0.38, 0.33, 0.28 fm-3 (2.2$\rho_{0}$,
1.9$\rho_{0}$, 1.6$\rho_{0}$) in the cases of $\Lambda 380$, $\Lambda 350$ and
$\Lambda 320$, respectively. Thin (thick) short-dashed curves are for the QHT
models with Q2 and Q3 (Q4). It should be noted that pressures in the
quarkyonic matter increase more rapidly with density than those in the QHT
matter. As discussed later, the rapid growth of pressure with density in the
range of $2\rho_{0}\sim 4\rho_{0}$ is an important feature of the quarkyonic
model. This rapid increase of pressure at onset of the quarkyonic phase
influences significantly on neutron-star $MR$ curves.

In Fig.3, pressures $P$ are drawn as a function of the energy density
$\varepsilon$, which are related closely to neutron-star $MR$ curves. The dot-
dashed (dotted) curve shows pressures in $\beta$-stable nucleonic (hadronic)
matter. Thin (thick) solid curves show pressures in quarkyonic matter for
$\Lambda 350$ and $\Lambda 320$ ($\Lambda 380$) with B1. The short-dashed
curve is for the QHT matter with Q4. Though the curves for Q4 and $\Lambda
380$ are rather similar to each other in comparison with the corresponding
curves in Fig.2, the former is still less steep than the latter in the region
of low energy density. As shown later, the EoSs for the QHT model Q4 and the
quarkyonic model $\Lambda 380$ lead to the neutron-star $MR$ curves consistent
with the criterion Eq.(1).

In Fig.4, sound velocities are drawn as a function of $\rho_{B}$. The dot-
dashed curve is sound velocities in $\beta$-stable nucleonic matter. Solid
curves are those in quarkyonic matter for $\Lambda 380$, $\Lambda 350$ and
$\Lambda 320$ with B1. There appear peak structures in the solid curves, being
related to rapid increasing of pressures in the range of $2\rho_{0}\sim
4\rho_{0}$. The dashed curve is sound velocities in the QHT matter with Q4 and
the dotted one is those in $\beta$-stable hadronic matter with hyperon mixing,
in which there appears no peak structure. The dashed curve becomes $c_{s}>c$
in high-density region. Also, the peak regions of solid curves become
$c_{s}>c$, if B2 is used instead of B1 for nuclear parts. In such regions of
$c_{s}>c$, sound velocities are approximated to be $c_{s}=c$.

It is interesting to notice that the peak structures in our quarkyonic-matter
results are somewhat similar to those for the QHC-matter EoS (QHC21) found in
KBH2022 . Our QHT-matter EoS gives no peak structure in sound velocities,
being different from both of them.

In the left panel of Fig.5, solid curves show pressures in quarkyonic matter
for $\Lambda 380$ in the cases of using B1 and B2 for nuclear interactions,
and short-dashed (dashed) curves are partial pressures of neutrons (quarks) in
respective cases. The dot-dashed curve is pressures in $\beta$-stable
nucleonic matter. Pressures in quarkyonic matter are found to be completely
dominated by neutron partial pressures. In order to reveal the reason why
neutron pressures in quarkyonic matter are far higher than those in
$\beta$-stable nucleonic matter, we show the neutron chemical potentials in
the cases of using B1 and B2 for nuclear interactions: In the right panel of
Fig.5, neutron chemical potentials $\mu_{n}$ are drawn as a function of
$\rho_{B}$. Lower and upper solid curves give neutron chemical potentials in
quarkyonic matter for $\Lambda 380$ in the cases of using B1 and B2,
respectively. The dot-dashed curve gives neutron chemical potential in
$\beta$-stable nucleonic matter. The neutron chemical potentials in quarkyonic
matter are far higher than those in the $\beta$-stable nucleonic matter, which
makes neutron pressures in the former far higher than those in the latter. The
reason of higher chemical potentials in the quarkyonic matter is because the
existence of free quarks inside the Fermi sea gives nucleons extra kinetic
energies by pushing them to higher momenta MR2019 .

###  III.2 $MR$ diagrams

We have the two types of hybrid EoSs, the QHT-matter EoS and the quarkyonic-
matter EoS. They are combined with the $\beta$-stable nucleonic-matter EoS
connected smoothly to the crust EoS Baym1 ; Baym2  in the low-density side.
The $MR$ relations of hybrid stars can be obtained by solving the TOV
equations with these hybrid EoSs.

Figure 6: (Color online) Star masses as a function of radius $R$. The dot-
dashed curves are by the $\beta$-stable nucleonic matter EoS. In the left
panel, thin (thick) solid curves are by the QHT-matter EoSs with Q2 and Q3
(Q4). The dotted curve is by the hadronic matter EoS including hyperons. In
the right panel, thin (thick) solid curves are by the quarkyonic-matter EoSs
for $\Lambda 350$ and $\Lambda 320$ ($\Lambda 380$) with B1. In both panels,
the horizontal lines indicates $R_{1.4M_{\odot}}=12.56^{+1.00}_{-1.07}$ km and
$R_{2.0M_{\odot}}=12.41^{+1.00}_{-1.10}$ km, and the rectangle indicates the
region of mass $M_{max}/M_{\odot}=2.21^{+0.31}_{-0.21}$ Legred2021 .

In Fig.6, star masses are given as a function of radius $R$. The dot-dashed
curves are obtained by the $\beta$-stable nucleonic matter EoS. In the left
panel, thin (thick) solid curves are obtained by the QHT-matter EoSs with Q2
and Q3 (Q4). The dotted curve is by the hadronic-matter EoS including
hyperons. In the cases of Q2, Q3 and Q4, the maximum masses are
$M_{max}/M_{\odot}$= 2.23, 2.30, 2.40, respectively, and the radii at
$2.0M_{\odot}$ are 11.8 km, 12.2 km, 12.5 km, respectively. In the right
panel, thin (thick) solid curves are obtained by the quarkyonic-matter EoSs
for $\Lambda 350$ and $\Lambda 320$ ($\Lambda 380$) with use B1 for nuclear
interactions. In the cases of $\Lambda 380$, $\Lambda 350$ and $\Lambda 320$,
the maximum masses are $M_{max}/M_{\odot}$= 2.64, 2.79, 2.76, respectively,
and the radii at $2.0M_{\odot}$ are 12.6 km, 13.1 km, 13.5 km, respectively.
In both panels, the horizontal lines indicates
$R_{1.4M_{\odot}}=12.56^{+1.00}_{-1.07}$ km and
$R_{2.0M_{\odot}}=12.41^{+1.00}_{-1.10}$ km, and the rectangle indicates the
region of mass $M_{max}/M_{\odot}=2.21^{+0.31}_{-0.21}$ Legred2021 . The thick
solid curve for Q4 in the left panel and that for $\Lambda 380$ in the right
panel are found to be consistent with the criterion Eq.(1), and the key
features of $R_{2M_{\odot}}\approx R_{1.4M_{\odot}}$ are found in these cases.

Then, it should be noted that the maximum mass $2.64M_{\odot}$ for $\Lambda
380$ is substantially larger than the value $2.40M_{\odot}$ for Q4. The reason
for such a difference between maximum masses can be understood by comparing
the $P(\rho_{B})$ curves in Fig.2, where the solid curve for $\Lambda 380$
increases more rapidly at onset of the quakyonic matter than the dashed curve
for Q4 at onset of quark matter. This means that the stiffness for former is
larger than that for the latter. In the case of QHT matter, it is not possible
to obtain such a rapid increasing of $P(\rho_{B})$ in the low-density region,
even if the $QQ$ repulsions are strengthened.

In the case of hadronic (nucleonic) matter, shown by the dotted (dot-dashed)
curve in the left panel, the maximum mass is 1.82$M_{\odot}$
(2.19$M_{\odot}$). The reduction of 0.37$M_{\odot}$ is due to the EoS
softening by hyperon (($\Lambda$ and $\Sigma^{-}$) mixing. This softening is
mainly caused by $\Sigma^{-}$ mixing: If only $\Lambda$ mixing is taken into
account, the maximum mass is obtained as 2.06$M_{\odot}$ being close to the
value of 2.19$M_{\odot}$ without hyperon mixing (dot-dashed curve). Thus,
massive stars with $M>2M_{\odot}$ cannot be obtained by the hadronic matter
EoSs with hyperon ($\Lambda$ and $\Sigma^{-}$) mixing YFYR14 ; YFYR16 ;
YTTFYR17 . On the other hand, the value of $R_{1.4M_{\odot}}$ is 12.4 (12.5)
km in the case of hadronic (nucleonic) matter, which means that the hyperon
mixing does not depend much on $R_{1.4M_{\odot}}$.

Figure 7: (Color online) Star masses as a function of central baryon density
$\rho_{Bc}$. The dot-dashed curves are by the $\beta$-stable nucleonic EoS.
The solid curve is by the quarkyonic-matter EoS for $\Lambda 380$ in the case
of using B1. The short-dashed curve is by the QHT-matter EoS for Q4.

In Fig.7, star masses are given as a function of central baryon density
$\rho_{Bc}$. The dot-dashed curves are by the $\beta$-stable nucleonic matter
EoS. The solid curve is obtained by the quarkyonic-matter EoS for $\Lambda
380$ with B1, and the dashed curve is by the QHT-matter EoS for Q4, where the
onset density in the former (latter) 0.39 (0.33) fm-3. Both of them are
consistent with Eq.(1), but the former mass curve for $\rho_{Bc}$ is
considerably above the latter one, as well as the corresponding $MR$ curves.

In Fig.8, star masses are given as a function of radius $R$. The solid curve
is obtained by the quarkyonic-matter EoS for $\Lambda 380$ with use of B1
($V_{nn}$) for nuclear interactions, given also in Fig.6. Dashed and short-
dashed curves are by the quarkyonic-matter EoSs for $\Lambda 380$ and $\Lambda
400$, respectively, in the case of using B2 ($V_{nn}$+$V_{nnn}$) instead of
B1. The difference between solid and dashed curves demonstrates the effect of
the three-neutron repulsion $V_{nnn}$, giving the larger maximum mass and
larger value of $R_{2.0M_{\odot}}$. The short-dashed curve for $\Lambda 400$
indicates that this effect of $V_{nnn}$ to increase mass and radius is
cancelled out by taking larger values of $\Lambda$.

Figure 8: (Color online) Star masses as a function of radius $R$. The dot-
dashed curves are by the $\beta$-stable nucleonic matter EoS. The solid curve
is for $\Lambda 380$ with B1. Dashed and short-dashed curves are by the
quarkyonic-matter EoSs for $\Lambda 380$ and $\Lambda 400$ with B2,
respectively. The horizontal dotted lines indicates
$R_{1.4M_{\odot}}=12.56^{+1.00}_{-1.07}$ km and
$R_{2.0M_{\odot}}=12.41^{+1.00}_{-1.10}$ km.

In Fig.9, star masses are given as a function of radius $R$. The solid curve
is obtained by the quarkyonic-matter EoS for $\Lambda 380$ with $\kappa=0.3$
in the case of using B1, given also in Fig.6. The dashed curve is obtained by
the approximation used in MR2019 , where the $QQ$ interactions are neglected
and the quark energy density Eq.(LABEL:edenq) is replaced by the kinetic
energy density. Then, the difference between short-dashed and dashed curves is
due to this approximation. The short-dashed curve is obtained by taking
$\kappa=0.4$ under this approximation. The similarity between solid and short-
dashed curves means that the deviation due to this approximation is canceled
out by adjusting the value of $\kappa$. In the same case of $\Lambda 380$ and
$\kappa=0.3$ with B1, the dotted curve is obtained by replacing the potential
energy density in Eq.(LABEL:edenn) to Eq.(22), being the approximated
treatment in MR2019 . This approximation to use Eq.(22) is found to reduce
masses and to increase radii.

Figure 9: (Color online) Star masses as a function of radius $R$. The dot-
dashed curves are by the $\beta$-stable nucleonic matter EoS. The solid curve
is obtained by the quarkyonic-matter EoS for $\Lambda 380$ with $\kappa=0.3$
in the case of using B1. The dashed (short-dashed) curve is for $\Lambda 380$
with $\kappa=0.3$ ($\kappa=0.4$) by the approximation to neglect potential
sectors in quark energy densities. The dotted curve is obtained by replacing
the potential energy density in Eq.(LABEL:edenn) to Eq.(22). The horizontal
lines indicates $R_{1.4M_{\odot}}=12.56^{+1.00}_{-1.07}$ km and
$R_{2.0M_{\odot}}=12.41^{+1.00}_{-1.10}$ km.

###  III.3 Discussion

In Legred2021 , they present the neutron-star properties such as maximum mass,
radius, tidal deformability, pressure and central density inferred from their
analysis, for which the median and $90\%$ highest-probability-density credible
regions are given. From Table II of Legred2021 , we choose the quantities in
the case of w/J0740+6620 Miller+ in order to compare with the corresponding
values obtained from our QHT-matter and the quarkyonic matter EoSs. In Table
1, tabulated are maximum masses $M_{max}$, pressures $p$ at $\rho_{0}$,
$2\rho_{0}$ and $6\rho_{0}$, radii $R$ and dimensionless tidal deformabilities
$\Lambda$ at $1.4M_{\odot}$ and $2.0M_{\odot}$, central densities $\rho_{c}$
at $1.4M_{\odot}$, $2.0M_{\odot}$ and $M_{max}$. Here, our results are for the
$\beta$-stable nucleonic matter EoS denoted as NUC, the QHT-matter EoS Q4 and
the quarkyonic matter EoS V380. These EoSs are adjusted so as to reproduce
$R_{1.4M_{\odot}}$ with an accuracy of a few hundred meters. Then, the key
feature of $R_{2M_{\odot}}\approx R_{1.4M_{\odot}}$ is found in the cases of
Q4 and V380 EoSs, contrastively to the case of the nucleonic EoS giving
$R_{2M_{\odot}}<R_{1.4M_{\odot}}$. The values of $R_{2.0M_{\odot}}$, central
densities and tidal deformabilities for Q4 and V380 EoSs are far closer to the
median values than those for nucleonic EoS, demonstrating the clear impacts of
quark phases in Q4 and V380 EoSs. The deviations from the median values in the
latter are considerably larger than those in the formers. Especially, the
values of $\Lambda_{1.4}$ and $\Lambda_{2.0}$ for the nucleonic EoS are noted
to be out of $90\%$ credible regions.

In the case of the quarkyonic matter EoS for V380, the values of $M_{max}$ and
$p(6\rho_{0})$ are found to be far larger than that for the nucleonic EoS. It
is interesting that such a large value of $M_{max}$ can be obtained
straightforwardly from the quarkyonic-matter EoS, considering the implication
of the large mass $(2.35\pm 0.17)M_{\odot}$ for PSR J0952-0607 Romani2022 .
The reason why a large value of $M_{max}$ is obtained n the case of the
quarkyonic matter EoS is because the pressure rises rapidly in the region of
$\rho_{B}\sim 2\rho_{0}$ as found in Fig.2. In the McLerran-Reddy model of the
quarkyonic matter, the resulting EoS is mainly controlled by the one parameter
$\Delta_{qyc}$ for Fermi-layer thickness. Then, it is difficult to reproduce
simultaneously $M_{max}=2.2M_{\odot}$ and $R_{2.0M_{\odot}}=12.4$ km.

Table 1: Maximum masses $M_{max}$, pressures $p$ at $\rho_{0}$, $2\rho_{0}$ and $6\rho_{0})$, radii $R$ and tidal deformabilities $\Lambda$ at $1.4M_{\odot}$ and $2.0M_{\odot}$, central densities $\rho_{c}$ at $1.4M_{\odot}$, $2.0M_{\odot}$ and $M_{max}$. Results for the $\beta$-stable nucleonic matter EoS denoted as NUC, the QHT-matter EoS Q4 and the quarkyonic matter EoS V380 are compared with the values taken from Legred2021 .  |  NUC |  Q4 |  $\Lambda 380$ | Ref.Legred2021
---|---|---|---|---
$M_{max}/M_{\odot}$ | 2.19 | 2.40 | 2.64 | $2.21^{+0.31}_{-0.21}$
$p(\rho_{0})$ ($10^{33}{\rm dyn/cm^{2}}$) | 5.27 | 5.27 | 5.27 | $4.30^{+3.37}_{-3.80}$
$p(2\rho_{0})$ ($10^{34}{\rm dyn/cm^{2}}$) | 2.76 | 5.09 | 4.42 | $4.38^{+2.46}_{-2.96}$
$p(6\rho_{0})$ ($10^{35}{\rm dyn/cm^{2}}$) | 6.94 | 12.0 | 22.6 | $7.41^{+5.87}_{-4.18}$
$R_{1.4M_{\odot}}$ (km) | 12.5 | 12.7 | 12.5 | $12.56^{+1.00}_{-1.07}$
$R_{2.0M_{\odot}}$ (km) | 11.8 | 12.5 | 12.6 | $12.41^{+1.00}_{-1.10}$
$R_{2.0M_{\odot}}-R_{1.4M_{\odot}}$ (km) | $-0.72$ | $-0.14$ | $+0.03$ | $-0.12^{+0.83}_{-0.85}$
$\Lambda_{1.4}$ | 779 | 525 | 473 | $507^{+234}_{-242}$
$\Lambda_{2.0}$ | 128 | 46 | 49 | $44^{+34}_{-30}$
$\rho_{c}(1.4M_{\odot})$ ($10^{14}{\rm g/cm^{3}}$) | 7.9 | 6.6 | 6.8 | $6.7^{+1.7}_{-1.3}$
$\rho_{c}(2.0M_{\odot})$ ($10^{14}{\rm g/cm^{3}}$) | 12. | 9.1 | 8.0 | $9.7^{+3.6}_{-3.1}$
$\rho_{c}(M_{max})$ ($10^{15}{\rm g/cm^{3}}$) | 1.8 | 1.6 | 1.3 | $1.5^{+0.3}_{-0.4}$

##  IV Conclusion

The observed masses and radii of neutron stars give constraints on the dense
matter EoSs and resulting $MR$ diagrams. In this sense, the observations of
massive stars over $2M_{\odot}$ and the NICER implication of
$R_{2M_{\odot}}\approx R_{1.4M_{\odot}}$ are critically important for
restricting neutron-star matter EoSs. In the case of hadronic matter, even if
the nucleonic matter EoS is constructed so as to be stiff enough to give the
maximum mass over $2M_{\odot}$, the hyperon mixing brings about a remarkable
softening of the EoS. The EoS-softening by hyperon mixing can be reduced, for
instance, by introducing many-body repulsions which work universally for every
kind of baryons. However, such a repulsive effect does not cancel out
completely the EoS softening by hyperon mixing: In the case of hadronic matter
EoS with hyperon mixing, it is difficult to obtain maximum masses over
$2M_{\odot}$. The most promising approach to solve this “hyperon puzzle” is to
assume the existence of quark phases in inner cores of neutron stars, namely
hybrid stars having quark matter in their cores.

When quark deconfinement phase transitions from a hadronic-matter EoS to a
sufficiently stiff quark-matter EoS are taken into account in the neutron-star
interiors, repulsive effects such as $QQ$ repulsions in quark phases are
needed in order to obtain sufficiently stiff EoSs resulting in massive hybrid
stars with masses over $2M_{\odot}$. In our QHT matter, it is possible to
reproduce maximum masses over $2M_{\odot}$ consistently with the NICER
implication, where the $QQ$ repulsion is taken to be strong enough and the
quark-hadron transition density is adjusted so as to be about $2\rho_{0}$ by
tuning of the density dependence of effective quark mass.

In the quarkyonic matter, the degrees of freedom inside the Fermi sea are
treated as quarks, and nucleons exist at the surface of the Fermi sea. The
existence of free quarks inside the Fermi sea gives nucleons extra kinetic
energy by pushing them to higher momenta. This mechanism of increasing
pressure is completely different from the above mechanism of EoS stiffening by
strong $QQ$ repulsions in the QHT matter. In calculations of $MR$ diagrams
with the quarkyonic-matter EoS, the critical quantity is the thickness
$\Delta_{qyc}$ of Fermi layer controlled by the parameters $\Lambda$ and
$\kappa$. With the reasonable choice of these parameters, the $MR$ curves of
quarkyonic hybrid stars are obtained so as to be consistent with the NICER
implication.

As well as $R_{2.0M_{\odot}}$, central densities and tidal deformabilities are
inferred from the analysis of the NICER data. The QHT-matter and quarkyonic
EoSs can be adjusted so as to reproduce these inferred quantities far closer
to the median values than those for nucleonic matter EoS, demonstrating the
clear impacts of quark phases in these cases..

Thus, the reasonable $MR$ curves of neutron stars can be derived from both
QHT-matter and quarkyonic-matter EoSs, having completely different mechanisms
to stiffen EoSs. However, when both EoSs are adjusted so as to be consistent
with the NICER implication, the maximum mass for the quakyonic-matter EoS is
considerably larger than that for the QHT-matter EoS.

## Acknowledgments

The authors would like to thank D. Blaschke for valuable comments and fruitful
discussions.

## References

  * (1) P.B. Demorest, T. Pennucci, S.M. Ransom, M.S.E. Roberts, and J.W. Hessels, Nature (London) 467, 1081 (2010).
  * (2) J. Antoniadis et al., Science 340, 6131 (2013).
  * (3) H.T. Cromartie et al., Nat. Astron. 4, 72 (2020).
  * (4) R.W. Romani et al., Astrophysical J. Lett. 934:L17 (2022).
  * (5) M.C. Miller, et al.. Astrophys. J. Lett. 918, L28 (2021), arXiv:2105.06979 [astro-ph.HE]
  * (6) T.E. Riley, et al.. Astrophys. J. Lett. 918, L27 (2021), arXiv:2105.06980 [astro-ph.HE]
  * (7) I. Legred, K. Chatziioannou, R. Essick, S. Han, and P. Landy, Phys. Rev. D 104. 063003 (2021),
  * (8) I. Bednarek, P. Haensel, J.L. Zdunik, M. Bejger, and R. Mańka, Astronomy & Astrophysics A157, 543 (2012).
  * (9) S. Weissenborn, D. Chatterjee, and J. Schaffner-Bielich, Nucl. Phys. A881, 62 (2012).
  * (10) M. Oertel, C.Providência, F. Gulminelli, and Ad.R. Raduta, J. Phys. G.42 075202 (2015).
  * (11) K.A. Maslov, E.E. Kolomeitsev, and D.N. Voskresensky, Phys. Lett. B748 369 (2015).
  * (12) S. Nishizaki, Y. Yamamoto, and T. Takatsuka, Prog. Theor. Phys. 105, 607 (2001); 108, 703 (2002).
  * (13) I. Vidan̂a,D. Logoteta, C. Providência, A. Polls, and I. Bombaci, Eur.Phys.Lett. 94, 1002 (2011).
  * (14) Y. Yamamoto, T. Furumoto. Yasutake, and Th.A Rijken, Phys Rev. C bf 90 045805 (2014).
  * (15) Y. Yamamoto, T. Furumoto, N. Yasutake, and Th.A Rijken, Eur. Phys. J. A52, 19 (2016).
  * (16) Y. Yamamoto, H. Togashi, T. Tamagawa. T. Furumoto, N. Yasutake, and Th.A. Rijken, Phys. Rev. C96, 065804 (2017).
  * (17) D. Lonardoni A. Lovato, S. Gandolfi, and F. Pederiva, Phys. Rev. Lett. 114, 092301 (2015).
  * (18) D. Logoteta, I. Vidan̂a, I. Bombaci, Eur. Phys. J. A55, 207 (2019).
  * (19) D. Cerstung, N. Kaiser. and W. Weise, Eur. Phys. J. A56, 175 (2020).
  * (20) A. Drago, A. Lavagno, G. Pagliara, and D. Pigato, Phys. Rev. C90, 065809 (2014).
  * (21) D.B. Kaplan, A.E. Nelson, Phys. Lett. B175, 57 (1986); 179, 409 (1986).
  * (22) G.E. Brown, C.-H. Lee, M. Rho, and V. Thorsson, Nucl. Phys. A567, 937 (1994).
  * (23) V. Thorsson, M. Prakash, and J.M. Lattimer, Nucl. Phys. A572, 693 (1994).
  * (24) C.-H. Lee, Phys. Rep. 275, 255 (1996).
  * (25) N.K. Glendenning and J. Schaffner-Bielich, Phys. Rev. Lett. 81, 4564 (1998).
  * (26) M.Baldo, G.F. Burgio, and H-J. Schulze, Superdense QCD Matter and Compact stars (NATO Science Series 2: Mathematics, Physics and Chemistry vol.197) (2006) ed D.Blaschke and D. Sedrakian.
  * (27) F. Özel, D. Psaltis, S. Ransom, P. Demorest, and M. Alford, Astrophys. J. Lett. 724, L199 (2010).
  * (28) I. Sagert, G. Pagliara, M, Hempel, and Schaffner-Bielich, Astrophys. J. Lett. 740, L14 (2011).
  * (29) T. Kähn, D. Blaschke, and D. Lastowiecki, Phys. Rev. D88, 085001 (2013).
  * (30) L. Bonanno and A. Sedrakian, Astrophys. J539, 416 (2012).
  * (31) R. Lastowiecki, D. Blaschke, H. Grigorian, and S. Typel, Acta Phys. Pol. B5, 535 (2012).
  * (32) M. Shahrbaf, D. Blaschke, A.G. Grunfeld, and H.R. Moshfegh, Phys. Rev. C101, 025807 (2020).
  * (33) M. Shahrbaf, D. Blaschke, and S. Khanmohamadi, J. Phys. G: Nucl. Part. Phys. 47, 115201 (2020).
  * (34) K. Otto, M. Oertel, and B-J Schaefer, Phys. Rev. D101, 103021 (2020); arXiv:1910.11929 (2020)
  * (35) T. Kojo, G. Baym, and T. Hatsuda, Astrophys. J. 934, 46 (2022), arXiv:2111.11919 (2022).
  * (36) Y. Yamamoto, N. Yasutake, and Th. A. Rijken, Phys. Rev. C105, 015804 (2022).
  * (37) A. Akmal, V.R. Pandharipande, and D.G. Ravenhall, Phys. Rev. C58, 1804 (1998).
  * (38) A.B. Kaidalov and K.A. Ter-Materosian, Nucl. Phys. B75, 471 (1974).
  * (39) M. M. Nagels, Th. A. Rijken, and Y. Yamamoto, Phys. Rev. C99, 044002 (2019).
  * (40) M. M. Nagels, Th. A. Rijken, and Y. Yamamoto, Phys. Rev. C99, 044003 (2019).
  * (41) M. M. Nagels, Th. A. Rijken, and Y. Yamamoto, Phys. Rev. C102, 054003 (2020).
  * (42) H. Togashi and M. Takano, Nucl. Phys. A902, 53 (2013).
  * (43) L. McLerran and H.D. Pisarski, Nucl. Phys. A796, 83 (2007).
  * (44) Y. Hidaka, L. McLerran, and H.D. Pisarski, Nucl. Phys. A808, 117 (2008).
  * (45) L. McLerra and S. Reddy, Phys. Rev. Lett. 122, 122701 (2019).
  * (46) S. Han, M.A.A. Mamun, S. Lalit, C. Constantinou, and M. Prakash, Phys. Rev. D100, 103022 (2019).
  * (47) D.C. Duarte, S. Hernande-Ortiz, and K.S. Jeong, Phps. Rev. C102, 025203 (2020); 065202 (2020).
  * (48) T. Zhao and J.M. Lattimer, Phys. Rev. D102, 023021 (2020).
  * (49) J. Margueron, H. Hansen, and P. Proust, and G.Chanfray, Phys. Rev. C104, 055803 (2021).
  * (50) G. Cao Phys. Rev. D105, 114020 (2022).
  * (51) N.K. Glendenning, Phys. Rev. D46, 1274 (1992).
  * (52) G. Baym, A. Bethe, and C. Pethick, Nucl. Phys. A175, 225 (1971).
  * (53) G. Baym, C.J. Pethick, and P. Sutherland, Astrophys. J.170, 299 (1971).