# Hybrid Star Properties with NJL and MFTQCD Model: A Bayesian Approach

Milena Albino<EMAIL_ADDRESS>Department of Physics, CFisUC,
University of Coimbra, P-3004 - 516 Coimbra, Portugal  Tuhin Malik<EMAIL_ADDRESS>Department of Physics, CFisUC, University of Coimbra, P-3004 - 516 Coimbra,
Portugal  Márcio Ferreira<EMAIL_ADDRESS>Department of Physics,
CFisUC, University of Coimbra, P-3004 - 516 Coimbra, Portugal  Constança
Providência<EMAIL_ADDRESS>Department of Physics, CFisUC, University of Coimbra,
P-3004 - 516 Coimbra, Portugal

###### Abstract

Context: The composition of the core of neutron stars is still under debate.
One possibility is that because of the high densities reached in their cores,
matter could be deconfined into quark matter.

Aims: In this work, we investigate the possible existence of hybrid stars,
using microscopic models to describe the different phases of matter. Within
the adopted microscopic models we aim at calculating both the properties of
neutron stars and the properties of matter. In particular, we want to probe
the influence of pQCD calculations and analyze the properties that identify a
transition to deconfined matter.

Methods: A Bayesian approach using a Markov-Chain Monte Carlo sampling process
is applied to generate 8 sets of equations of state. A Maxwell construction is
adopted to describe the deconfinement transition. For the hadronic phase, we
consider a stiff and a soft EOS obtained from the Relativistic Mean Field
(RMF) model with non-linear meson terms. For the quark phase, we use two
different models: the Nambu-Jona-Lasinio (NJL) model with multiquark
interactions and the Mean Field Theory of QCD, a model similar to the MIT bag
model with a vector term. Bayesian inference was applied to determine the
model parameters that satisfy the X-ray observations from NICER and have phase
transition at densities between 0.15 and 0.40 fm-3. We have also applied
restrictions from the pQCD calculations to half of the sets.

Results: Hybrid stars are compatible with current observational data. The
restrictions of pQCD reduce the value of the maximum mass. However, even when
applying this restriction, the models were able to reach values of
$M_{\text{max}}=2.1-2.3M_{\odot}$. The conformal limit was still not attained
at the center of the most massive stars. The vector interactions are essential
to describe hybrid stars with a mass above 2$M_{\odot}$. The multiquark
interactions introduced may affect the limits of some quantities considered as
indicators of the presence of a deconfined phase. It is possible to find a set
of EOS, that predict that inside NS the renormalized matter trace anomaly is
always positive.

##  I Introduction

Neutron stars (NSs) represent some of the densest forms of matter observable
in the universe, providing unique natural laboratories for studying the
extreme states of quantum chromodynamics (QCD) [1, 2]. The cores of NSs, in
particular, may host a myriad of exotic phases of matter, offering a window
into the behavior of QCD under such extreme conditions [3, 4]. Theoretical
predictions based on perturbative QCD (pQCD) suggest the presence of
deconfined quark matter at extraordinarily high baryon densities,
approximately forty times greater than the density of nuclear matter at
saturation $\rho_{0}$ [5]. The most challenging region to study theoretically
is, however, at intermediate densities i.e. few times nuclear matter
saturation density which is actually relevant for the matter in the core of
NSs [3].

At these intermediate densities, fundamental calculations like Lattice QCD
face significant computational obstacles, particularly the well-known sign
problem in simulations at finite baryon densities. The current ability of
chiral effective field theory ($\chi$EFT) to accurately compute the EOS at
relatively low baryon densities, below twice saturation density [6, 7] does
not cover the range of densities that is hypothesized in neutron star (NS)
cores. In contrast, various effective models that suggest a range of novel and
exotic phases of quark matter in the intermediate density have been studied
over the past decades. These include pion superfluidity [8, 9, 10], different
types of color superconductivity such as two-flavor color superconductivity
[11, 12, 13, 14] and the color flavor locked (CFL) phase [15], along with more
intricate configurations like the Larkin-Ovchinkov-Fulde-Ferrel (LOFF) [16,
17] and crystalline superconductivity phases. To support the possibility of
quark matter cores within neutron stars, recent model-independent studies
combining astrophysical observations and theoretical calculations have shown
that the interior properties of the most massive neutron stars align with the
characteristics of a deconfined quark phase, suggesting the existence of
sizable quark-matter cores in stars approaching or exceeding two solar masses
[18, 19, 20]. However, no information on the composition was considered in the
intermediate density regime. The behavior of some quantities as the speed of
sound, the polytropic index [18], the renormalized trace anomaly
$\Delta=1/3-P/\epsilon$ [21] and the measure of conformality
$d_{c}\equiv\sqrt{\Delta^{2}+\Delta^{\prime 2}}$, where $\Delta^{\prime}$ is
the logarithmic derivative with respect to the energy density [22], were
analyzed and proposed as possibly indicating the existence of a quark core in
NSs. These analyzes were carried out considering agnostic model independent
descriptions of the EOS.

One of the aims of this work is to describe hybrid stars considering
physically justified hadron and quark models in order to study the possible
existence of a quark core inside a NS. The quantities that have been proposed
to identify the presence of a quark core within agnostic models will be
studied considering these microscopic hybrid EOS, and confronted with the
values proposed as indicative of the presence of deconfined matter.

Since the formulation of QCD, the possible existence of quarks inside neutron
stars has been questioned [3]. One of the first models of hybrid stars [23]
predicted that the quark-hadron transition would only happen around
$10\rho_{0}$, making the existence of hybrid stars very unlikely. However, as
new models emerged, the results pointed to a phase transition at lower
densities, in which it was possible to have quarks inside neutron stars. Since
then, various studies have been carried out, considering different hadronic
models [24, 25, 26, 6], hybrid models [13, 27, 28, 29, 30, 31, 32, 33, 34, 35,
36, 37, 38, 39, 40, 41] or even model-independent studies [18, 19, 21, 22].
Despite numerous studies and the rapid development of pulsar detection
techniques, it is still not possible to determine with confidence the
existence of quarks in the NS. In this work, we apply a Bayesian inference
approach to generate numerous hybrid equations of state (EOS) that satisfy the
observational data, in order to explore the possible existence of hybrid
stars. We adopt the two-model approach, in which the Relativistic Mean Field
(RMF) model, as in [25], describes the hadron phase, and the quiral symmetric
SU(3) Nambu-Jona-Lasinio (NJL) model [42, 43, 44] with four-quark and eight-
quark interaction terms [45, 29] describes the quark phase. We also obtain
hybrid EOS sets using the Mean Field Approximation of QCD (MFTQCD). The MFTQCD
model is derived from the QCD Lagrangian, considering a gluon field
decomposition in soft and hard momentum components. A mean field approximation
is considered for the hard gluons and its contribution stiffens the EOS
allowing for high pressures at high energy densities. The soft gluon fields
originate condensates that soften the EOS and give rise to a bag like term
[46].

In 1961, before the formulation of QCD, Yochiro Nambu and Giovanni Jona-
Lasinio proposed the Nambu-Jona-Lasinio (NJL) model [47, 48]. In the original
idea, the Dirac particle mass was assumed to be produced by the nucleon
interaction, analogous to the superconductivity in BCS theory. Later, NJL
model was adapted to have quarks and gluons as degrees of freedom. It became
an effective theory of QCD, being a good approximation in the low energy
limit. In the NJL model, gluons are omitted and the quark interaction is
assumed to be point-like. Although color confinement property can not be
described by NJL, it reproduces all the QCD global symmetries and the
spontaneus chiral symmetry breaking, see [49] for a review. The NJL model is
suitable in the context of NS and has been widely used to describe the
deconfined phase in hybrid stars [44, 50, 51, 52, 53, 54, 13, 55, 29, 45].

Some of the studies discussing hybrid stars also apply a Bayesian inference to
constrain the model parameters, either using a grid approach [31, 38, 40, 41]
or a Markov-Chain Monte Carlo sampling process [39]. Different models are
used, a metamodeling description for the hadronic [38, 39], or a RMF model
[31] together with a quark phase described by the NJL model [31, 38] or a
quark metamodel [39], in all cases describing the deconfinement phase
transition by a Maxwell construction, or a parametrized smooth crossover
between the hadronic phase described by a RMF EOS and a quark phase described
within a constituent quark model [41].

In a previous work [25], Bayesian inference was used to explore the NSs EOS
with a RMF model. However, it was assumed that NS were made entirely by hadron
matter. In this work, we explore the possibility of the existence of quark
matter inside NS using microscopic models. For the hadron phase, we choose two
equations of state from [25] with extreme properties, in particular, a soft
and a stiff EOS. The quark phase is described by the Nambu-Jona-Lasinio (NJL)
model and the MFTQCD model. Bayesian inference is used to obtain large samples
of EOS that satisfy the imposed constraints, which include the NICER x-ray
observations from PSR J0030+0451 [56, 57] and PSR J0740+6620 [58, 59] and a
constraint in the phase transition density, which we consider should occur
above saturation density. In [60] it was found that pQCD calculations can have
an influence at intermediate densities, if thermodynamic relations and
causality are taken into account. Therefore, we also impose these pQCD
constraints in some of the generated sets to analyze their effect on the
results for NSs.

This article is structured as follows: Section II gives a brief overview of
the NJL, MFTQCD and RMF models. The Bayesian inference is described in Section
III. The results are presented and discussed in Section IV while the
conclusions are drawn in Section V.

##  II Model

In this section, we present the quark and the hadron models separately. To
build hybrid equations, we assumed a first order quark-hadron phase transition
and applied the Maxwell construction.

###  II.1 Quark phase

We obtained sets considering the Nambu-Jona-Lasinio (NJL) model (Section
II.1.1) and sets considering the Mean Field Theory of QCD (Section II.1.2). A
brief review of both models is presented here.

####  II.1.1 Nambu-Jona-Lasinio (NJL)

In this work, we considered the SU(3) Nambu-Jona-Lasinio (NJL) Lagrangian
being as

$\displaystyle\mathcal{L}_{\text{NJL}}$
$\displaystyle=\bar{\psi}\left(i\not{\partial}-m+\mu\gamma^{0}\right)\psi$
$\displaystyle\quad+\frac{G}{2}\left[\left(\bar{\psi}\lambda_{a}\psi\right)^{2}+\left(\bar{\psi}i\gamma^{5}\lambda_{a}\psi\right)^{2}\right]$
$\displaystyle\quad+\mathcal{L}_{\text{'t Hooft}}+\mathcal{L}_{\text{I}},$ (1)

which $\mathcal{L}_{\text{'t Hooft}}$ is known as the ’t Hooft term, to
simulate the U(1)A symmetry breaking

$\mathcal{L}_{\text{'t
Hooft}}=8\kappa\left[\det\left(\bar{\psi}P_{R}\psi\right)+\det\left(\bar{\psi}P_{L}\psi\right)\right],$
(2)

where the determinant is applied in the space of flavors. However, these terms
are not enough to reach the $2M_{\odot}$. For this reason, we add 4- and
8-quark interaction terms:

$\displaystyle\mathcal{L}_{\text{I}}$
$\displaystyle=-G_{\omega}\left[\left(\bar{\psi}\gamma^{\mu}\lambda_{0}\psi\right)^{2}+\left(\bar{\psi}\gamma^{\mu}\gamma_{5}\lambda_{0}\psi\right)^{2}\right]$
$\displaystyle\quad-
G_{\rho}\sum_{a=1}^{8}\left[\left(\bar{\psi}\gamma^{\mu}\lambda^{a}\psi\right)^{2}+\left(\bar{\psi}\gamma^{\mu}\gamma_{5}\lambda^{a}\psi\right)^{2}\right]$
$\displaystyle\quad-
G_{\omega\omega}\left[\left(\bar{\psi}\gamma^{\mu}\lambda_{0}\psi\right)^{2}+\left(\bar{\psi}\gamma^{\mu}\gamma_{5}\lambda_{0}\psi\right)^{2}\right]^{2}$
$\displaystyle\quad-
G_{\sigma\omega}\sum_{a=0}^{8}\left[\left(\bar{\psi}\lambda_{a}\psi\right)^{2}+\left(\bar{\psi}i\gamma^{5}\lambda_{a}\psi\right)^{2}\right]$
$\displaystyle\qquad\times\left[\left(\bar{\psi}\gamma^{\mu}\lambda_{0}\psi\right)^{2}+\left(\bar{\psi}\gamma^{\mu}\gamma_{5}\lambda_{0}\psi\right)^{2}\right]$
$\displaystyle\quad-
G_{\rho\omega}\sum_{a=1}^{8}\left[\left(\bar{\psi}\gamma^{\mu}\lambda_{0}\psi\right)^{2}+\left(\bar{\psi}\gamma^{\mu}\gamma_{5}\lambda_{0}\psi\right)^{2}\right]$
$\displaystyle\qquad\times\left[\left(\bar{\psi}\gamma^{\mu}\lambda_{a}\psi\right)^{2}+\left(\bar{\psi}\gamma^{\mu}\gamma_{5}\lambda_{a}\psi\right)^{2}\right].$
(3)

Most of the interaction terms effects were studied in [55, 45, 29]. As the NJL
model is not renormalizable, we implemented the 3-momentum cutoff scheme.
After the mean field approximation, the grand canonical potential at $T=0$ can
be written as

$\displaystyle\Omega$
$\displaystyle=\Omega_{0}+G\left(\sigma_{u}^{2}+\sigma_{d}^{2}+\sigma_{s}^{2}\right)+4\kappa\sigma_{u}\sigma_{d}\sigma_{s}$
$\displaystyle\quad-\frac{2}{3}G_{\omega}(\rho_{u}+\rho_{d}+\rho_{s})^{2}$
$\displaystyle\quad-\frac{4}{3}G_{\rho}(\rho_{u}^{2}+\rho_{s}^{2}+\rho_{s}^{2}-\rho_{u}\rho_{d}-\rho_{u}\rho_{s}-\rho_{d}\rho_{s})$
$\displaystyle\quad-\frac{4}{3}G_{\omega\omega}(\rho_{u}+\rho_{d}+\rho_{s})^{4}$
$\displaystyle\quad-4G_{\sigma\omega}\left(\sigma_{u}^{2}+\sigma_{d}^{2}+\sigma_{s}^{2}\right)\left(\rho_{u}+\rho_{d}+\rho_{s}\right)^{2}$
$\displaystyle\quad-\frac{8}{3}G_{\rho\omega}\left(\rho_{u}+\rho_{d}+\rho_{s}\right)^{2}$
$\displaystyle\qquad\times\left(\rho_{u}^{2}+\rho_{d}^{2}+\rho_{s}^{2}-\rho_{u}\rho_{d}-\rho_{u}\rho_{s}-\rho_{d}\rho_{s}\right)$
$\displaystyle\quad-\frac{3}{\pi^{2}}\sum_{f=u,d,s}\int_{k_{F_{f}}}^{\Lambda}dpp^{2}E_{f}-\frac{1}{\pi^{2}}\sum_{f=u,d,s}\tilde{\mu}_{f}k_{F_{f}},$
(4)

where $k_{F_{f}}=\sqrt{\tilde{\mu}_{f}^{2}-\tilde{m}_{f}^{2}}$ is the Fermi
momentum and $\Omega_{0}$ is calculated to vanishs the potential in the
vacuum. The effective mass ($\tilde{m}$) and chemical potential
($\tilde{\mu}$) expressions are given by the gap equations

$\displaystyle\tilde{m}_{i}$
$\displaystyle=m_{i}+-2G\sigma_{i}-2\kappa\sigma_{j}\sigma_{k}$
$\displaystyle\quad+\frac{8}{3}G_{\sigma\omega}(\sigma_{i}^{2}+\sigma_{j}^{2}+\sigma_{k}^{2})\sigma_{i}$
(5) $\displaystyle\tilde{\mu}_{i}$
$\displaystyle=\mu_{i}-\frac{4}{3}G_{\omega}(\rho_{i}+\rho_{j}+\rho_{k})-\frac{4}{3}G_{\rho}(2\rho_{i}-\rho_{j}-\rho_{k})$
$\displaystyle\quad-\frac{16}{9}G_{\omega\omega}(\rho_{i}+\rho_{j}+\rho_{k})^{3}$
$\displaystyle\quad-\frac{8}{3}G_{\sigma\omega}(\sigma_{i}^{2}+\sigma_{j}^{2}+\sigma_{k}^{2})(\rho_{i}+\rho_{j}+\rho_{k})$
$\displaystyle\quad-\frac{8}{9}G_{\rho\omega}(\rho_{i}+\rho_{j}+\rho_{k})$
$\displaystyle\quad\times(4\rho_{i}^{2}+\rho_{j}^{2}+\rho_{k}^{2}-\rho_{i}\rho_{j}-\rho_{i}\rho_{k}-4\rho_{j}\rho_{k}),$
(6)

for ${i\neq j\neq k}\in{u,d,s}$. Imposing thermodynamic consistence [61],
i.e., $d\Omega/d\sigma_{f}=0$ and $d\Omega/d\rho_{f}=0$, we obtain expressions
for the quark condensate ($\sigma_{i}$) and density ($\rho_{i}$). At zero
temperature, we can write them as

$\displaystyle\sigma_{f}$
$\displaystyle=-\frac{3}{\pi^{2}}\int_{k_{F_{f}}}^{\Lambda}dpp^{2}\frac{\tilde{m}_{f}}{E_{f}},$
(7) $\displaystyle\rho_{f}$ $\displaystyle=\frac{1}{\pi^{2}}k_{F_{f}}^{3}.$
(8)

The parameters are set in order to reproduce the experimental values of mass
and decay constant from $\pi^{\pm}$, $K^{\pm}$, $\eta$ and $\eta^{\prime}$
mesons obtained in [62]. Their values are given in Table 1. In addition, we
also use the electron mass being $m_{e}=0.510\text{ MeV}$.

$\Lambda$ (MeV) |  $m_{u,d}$ (MeV) |  $m_{s}$ (MeV) | $G\Lambda^{2}$ | $\kappa\Lambda^{5}$
---|---|---|---|---
623.58 | 5.70 | 136.60 | 3.34 | -13.67
Table 1: Fixed parameters of NJL model. $\Lambda$ is the cutoff, $m_{u,d,s}$
are the quarks up, down and strange masses, $G$ and $\kappa$ are coupling
constants. The values of these parameters are obtained in such a way that they
reproduce the experimental masses and decay constants of mesons [62].

As the NJL is normalized in the vacuum, we add a constant $B$ in the pressure
equation, $P\to P+B$, similar to the MIT bag model. This term was also
included in [54, 13], in the first study to allow the restoration of chiral
symmetry and the deconfinement phase transition to occur at the same baryon
density, and in the second to control the deconfinement phase transition. In
the present study we include the extra $B$ term with the same proposal as in
[13]. In this way, we have 6 parameters to be determined by Bayesian analysis:
the adimensional coupling constants (defined as $\xi_{i}=2G_{i}/G$ and
$\xi_{ii}=16G_{ii}/G^{4}$) $\xi_{\omega}$, $\xi_{\rho}$, $\xi_{\omega\omega}$,
$\xi_{\omega\sigma}$, $\xi_{\omega\rho}$ and the constant $B$.

####  II.1.2 Mean Field Theory of QCD (MFTQCD)

In this model, we start with the QCD lagrangian and we assume we can decompose
the gluon field in low and high momentum components (soft and hard gluons)
[46, 63, 64]

$\tilde{G}^{a\mu}(k)=\tilde{A}^{a\mu}(k)+\tilde{\alpha}^{a\mu}(k),$ (9)

where $\tilde{G}$ is the gluon field in the momentum space and we defined $A$
and $\alpha$ as the soft and hard gluon, respectively. We approximate the
momentum of the soft gluons to be zero. This implies that their field is
constant. We replace the fields with their expected value in the vacuum [65,
66]

$\displaystyle\left<A^{a\mu}A^{b\nu}\right>$
$\displaystyle=-\frac{\delta^{ab}}{8}\frac{g^{\mu\nu}}{4}\mu_{0}^{2},$ (10)
$\displaystyle\left<A^{a\mu}A^{b\nu}A^{c\rho}A^{d\eta}\right>$
$\displaystyle=\frac{\phi_{0}^{4}}{(32)(34)}\left[g_{\mu\nu}g^{\rho\eta}\delta^{ab}\delta^{cd}\right.$
$\displaystyle\quad\left.+g_{\mu}^{\rho}g_{\nu}^{\eta}\delta^{ac}\delta^{bd}+g_{\mu}^{\eta}g_{\nu}^{\rho}\delta^{ad}\delta^{bc}\right],$
(11)

where $\phi_{0}$ and $\mu_{0}$ are the energy scale that will be determined.
Next, we consider the hard gluons as classical field, due to its large
occupation numbers at all energy levels [67, 68]

$\alpha_{\mu}^{a}\to\left<\alpha_{\mu}^{a}\right>=\alpha_{0}^{a}\delta_{\mu
0},$ (12)

where $\alpha_{0}^{a}$ is constant. After some calculation, we obtain the
MFTQCD lagrangian

$\displaystyle\mathcal{L}_{\text{MFTQCD}}$
$\displaystyle=-B+\frac{m_{G}^{2}}{2}\alpha_{0}^{a}\alpha_{0}^{a}$
$\displaystyle\quad+\sum_{q=1}^{N_{f}}\bar{\psi}_{i}^{q}\left(i\delta_{ij}\gamma^{\mu}\partial_{\mu}+g\gamma^{0}T_{ij}^{a}\alpha_{0}^{1}-\delta_{ij}m_{q}\right)\psi_{j}^{q},$
(13)

where we defined

$\displaystyle m_{G}^{2}$ $\displaystyle=\frac{9}{32}g^{2}\mu_{0}^{2},$ (14)
$\displaystyle B$
$\displaystyle=\frac{9}{4(34)}g^{2}\phi_{0}^{4}=\left<\frac{1}{4}F^{a\mu\nu}F^{a}_{\mu\nu}\right>.$
(15)

The 2-dimensional condensate contributes to the mass of the hard gluon
($m_{G}$), while the 4-dimensional one acts as a constant analogous to the MIT
bag model. After obtaining the equations of motion and using the energy-
momentum tensor, we obtain the equations of state below:

$\displaystyle P$ $\displaystyle=\frac{27}{2}\xi^{2}\rho_{B}^{2}-B+P_{F},$
(16) $\displaystyle\epsilon$
$\displaystyle=\frac{27}{2}\xi^{2}\rho_{B}^{2}+B+\epsilon_{F},$ (17)

where $P_{F}$ and $\epsilon_{F}$ are the pressure and energy density of a non-
interacting Fermi gas of quarks and electrons and $\xi=g/m_{G}$. The deduction
of the EOS can be found in [46]. We considered the masses values as
$m_{u}=5\text{ MeV}$, $m_{d}=7\text{ MeV}$ and $m_{s}=100\text{ MeV}$. $\xi$
and $B$ are determined from the Bayesian analysis.

###  II.2 Hadron phase

In [25], equations of state of pure nuclear matter were obtained from a
relativistic mean field model with non-linear meson terms, both self-
interactions and mixed terms. In this model, the nuclear interactions are
described by the exchange of the $\sigma$ (scalar-isoscalar), $\omega$
(vector-isoscalar) and $\rho$ (vector-isovector) mesons. The Lagrangian is
given by

$\mathcal{L}=\mathcal{L}_{N}+\mathcal{L}_{M}+\mathcal{L}_{NL},$ (18)

where

$\displaystyle\mathcal{L}_{N}$
$\displaystyle=\bar{\Psi}\left[\gamma^{\mu}\left(i\partial_{\mu}-g_{\omega}\omega_{\mu}-g_{\rho}\boldsymbol{t}\cdot\boldsymbol{\rho}_{\mu}\right)\right.$
$\displaystyle\quad\left.-(m-g_{\sigma}\phi)\right]\Psi,$ (19)
$\displaystyle\mathcal{L}_{M}$
$\displaystyle=\frac{1}{2}\left[\partial_{\mu}\phi\partial^{\mu}\phi-
m_{\sigma}^{2}\phi^{2}\right]$
$\displaystyle\quad-\frac{1}{4}F_{\mu\nu}^{(\omega)}F^{(\omega)\mu\nu}+\frac{1}{2}m_{\omega}^{2}\omega_{\mu}\omega^{\mu}$
$\displaystyle\quad-\frac{1}{4}\boldsymbol{F}_{\mu\nu}^{(\rho)}\cdot\boldsymbol{F}^{(\rho)\mu\nu}+\frac{1}{2}m_{\rho}^{2}\boldsymbol{\rho}_{\mu}\cdot\boldsymbol{\rho}^{\mu},$
(20) $\displaystyle\mathcal{L}_{NL}$
$\displaystyle=-\frac{1}{3}bmg_{\sigma}^{3}(\sigma)^{3}-\frac{1}{4}cg_{\sigma}^{4}(\sigma)^{4}+\frac{\xi}{4!}g_{\omega}^{4}(\omega_{\mu}\omega^{\mu})^{2}$
$\displaystyle\quad+\Lambda_{\omega}g_{\rho}^{2}\boldsymbol{\rho}_{\mu}\cdot\boldsymbol{\rho}^{\mu}g_{\omega}^{2}\omega_{\mu}\omega^{\mu},$
(21)

where the Dirac spinor $\Psi$ represents the nucleon doublet (neutron and
proton) with a bare mass $m$; $g_{\sigma}$ ($m_{\sigma}$), $g_{\omega}$
($m_{\omega}$) and $g_{\rho}$ ($m_{\rho}$) are the coupling constants (masses)
of the $\sigma$, $\omega$ and $\rho$ mesons, respectively; $b$, $c$, $\xi$ and
$\Lambda_{\omega}$ are the couplings that determine the strength of the non-
linear terms interaction. All coupling constants were obtained through a
Bayesian inference approach. This method was applied so that the equations
should satisfy the chiral EFT constraints for pure neutron matter, properties
of nuclear matter and neutron star observational data. In this work, we will
use two equations of state (EOS) from [25]: the BMPF 220, a soft EOS, and the
BMPF 260, a stiff EOS. Nuclear matter and NSs properties of these equations
are shown in Table 2. With this choice we analyze the extreme hadronic cases,
and consider they will define for a given mass the low and high radius
extremes.

Quantity | Units | BMPF220 | BMPF260
---|---|---|---
$\rho_{0}$ | fm-3 | 0.146 | 0.149
$m^{\ast}$ | … | 0.74 | 0.61
$\epsilon_{0}$ | MeV | -15.93 | -16.02
$K_{0}$ | MeV | 288 | 244
$J_{\text{sym},0}$ | MeV | 31.74 | 29.06
$L_{\text{sym},0}$ | MeV | 45 | 56
$K_{\text{sym},0}$ | MeV | -159 | 2
$M_{\text{max}}$ | M⊙ | 2.200 | 2.592
$R_{\text{max}}$ | km | 11.35 | 12.64
$R_{1.4}$ | km | 12.94 | 13.51
$\Lambda_{1.4}$ |  | 548 | 828
Table 2: Nuclear matter properties and NSs properties of the hadronic EOS.

##  III Bayesian Approach

The Bayesian method is a robust technique that enables us to sample an
ensemble of EOS, all of which meet the specified constraints. To apply this,
we have to define the prior, likelihood, and sampler.

Prior:- In Bayesian inference, a prior distribution represents our initial
assumptions about a parameter before any data is observed. This prior is then
revised with new data through the likelihood of the observed data to create
the posterior distribution. In our calculation, we have employed uniform
priors for all the parameters of the model. The minimum and maximum values of
the uniform distribution are shown in Table 3.

NJL |  | MFTQCD
---|---|---
Parameters | Units | min | max |  | Parameters | Units | min | max
$\xi_{\omega}$ | … | 0 | 0.6 |  | $\xi$ | MeV-1 | 0 | 0.0018
$\xi_{\rho}$ | … | 0 | 1 |  | $B$ | MeV.fm-3 | 50 | 180
$\xi_{\omega\omega}$ | … | 0 | 40 |  |  |  |  |
$\xi_{\sigma\omega}$ | … | 0 | 8 |  |  |  |  |
$\xi_{\rho\omega}$ | … | 0 | 50 |  |  |  |  |
$B$ | MeV.fm-3 | 0 | 50 |  |  |  |  |
Table 3: The lowest and highest values for the uniform distribution prior used
for NJL and MFTQCD in this study.

Likelihood:- In Bayesian inference, the likelihood quantifies the probability
of the observed data for particular parameter values, acting as a link between
prior beliefs and posterior knowledge.

  1. (i)

X-ray observation(NICER): X-ray observations give the mass and radius
measurements of NS. Therefore, the corresponding evidence takes the following
form,

$\displaystyle P(d_{\rm
X-ray}|\mathrm{EOS})=\int^{M_{\mathrm{max}}}_{M_{\mathrm{min}}}dmP(m|\mathrm{EOS})$
$\displaystyle\times P(d_{\rm X-ray}|m,R(m,\mathrm{EOS}))$
$\displaystyle=\mathcal{L}^{\rm NICER}$ (22)

$P(m|\rm{EOS})=\left\\{\begin{matrix}\frac{1}{M_{\mathrm{max}}-M_{\mathrm{min}}}&\text{
if }M_{\mathrm{min}}\leq m\leq M_{\mathrm{max}}\ ,\\\
0&\text{otherwise.}&\end{matrix}\right.$ (23)

Here, $M_{\rm min}$ is 1 $M_{\odot}$, and $M_{\rm max}$ represents the maximum
mass of a NS for the given EOS [69].

  2. (ii)

Phase transition: We constraint the phase transition value through a super-
Gaussian centered at $\rho_{\text{trans}}=0.275\text{ fm}^{-3}$, with a
standard deviation of $\sigma=0.08\text{ fm}^{-3}$ and $p=5$. Therefore, the
likelihood is written as

$\mathcal{L}^{\rm
PhT}=\text{exp}\left[-\left(\frac{(\rho_{\text{trans}}-\rho_{B})^{2}}{2\sigma^{2}}\right)^{p}\right].$
(24)

The density range was chosen such that the occurrence of the transition is
possible in medium mass stars with $M\geq 1.4M_{\odot}$. If the transition
occurs at a too large density the quark core is too small or does not exist
inside a NS. For the lower limit, saturation density was considered. Since we
are dealing with extremely asymmetric matter, with a proton fraction below 10%
it is reasonable to suppose that a transition to quark matter could be favored
energetically.

  3. (iii)

pQCD EOS: In [60], the authors found that for the neutron stars’ EOS to be in
agreement with the pQCD calculations, they must be within a limited region in
$P\times\epsilon\times\rho_{B}$ space. To apply this restriction, we have
calculated the area in $P\times\epsilon$ allowed by pQCD at
$\rho_{B}=7\rho_{0}$ for the renormalization scale $X=4$. Thus, we define the
likelihood as

$\mathcal{L}^{\rm
pQCD}=\begin{cases}1,(P(7\rho_{0}),\epsilon(7\rho_{0}))\text{ inside pQCD
region};\\\ 0,(P(7\rho_{0}),\epsilon(7\rho_{0}))\text{ {not}
inside},\end{cases}$ (25)

where $(P(7\rho_{0}),\epsilon(7\rho_{0}))$ is the pressure and energy density
at $7\rho_{0}$.

The final likelihood is then given by

$\mathcal{L}=\mathcal{L}^{\rm NICERI}\mathcal{L}^{\rm NICERII}\mathcal{L}^{\rm
PhT}\mathcal{L}^{\rm pQCD}\,.$ (26)

NICER I and NICER II correspond to the mass-radius measurements of PSR
J0030+0451 and PSR J0740+6620, respectively.

Sampler:- We utilized the nested sampling method, specifically the PyMultinest
[70, 71] sampler, as part of the Bayesian Inference Library (BILBY) [72].

##  IV Results

In this section we present and discuss the properties of hybrid stars,
considering different models to describe the total EOS. For the confined
phase, we choose two representative hadronic EOS, one soft and one hard, but
still compatible with the constraints imposed on the EOS. This choice allows
us to estimate the effect of the hadronic EOS on the total EOS. For the
deconfined phase, we consider the two different models, NJL and MFTQCD,
introduced in Sections II.1.1 and II.1.2, respectively. This choice allows us
to estimate the effects of the choice of a particular quark model. In order to
understand the effect of imposing the thermodynamic and casuality constraints
from the pQCD EOS, we have also performed the inference analysis with and
without this constraint in the following. Note that the NJL model cannot be
extended up to a density as large as the one where the pQCD EOS was
calculated, but following [60] we can estimate the effects of this constraint
on the hybrid EOS.

Figure 1: Pressure versus energy density for the soft EOS (left panel) and
stiff EOS (right panel). Gray plots show the results from [18].

Taking into account the above choices, we have generated, within a Bayesian
inference approach for the quark model parameters, eight different EOS sets
with $\sim 11000$ (NJL) and $\sim 2000$ (MFTQCD) EOS each combining the
following options: i) two hadron EOS (soft and hard), ii) two quark models
(NJL and MFTQCD), iii) with or without imposing the pQCD constraint. These
eight sets are discussed next.

In Fig. 1, the EOS distributions, i.e. the pressure versus the energy density,
are shown at a 90% confidence level (CL) for the two hadronic models, the soft
(stiff) on the left (right) panel. In both panels the hadronic EOS is
represented by a dotted line. The sets including (not including) the pQCD
constraint are represented by colored homogeneous bands (patterns) and the
sets with the NJL (MFTQCD) model are represented by the colors blue/dark blue
(orange/dark orange). The gray band represents the results from [18]. This
band spans the entire phase space generated by an agnostic description of the
EOS and imposing constraints on pure neutron matter from $\chi$EFT
calculations, the pQCD constraints at very high densities and NS observables.
Some conclusions are in order: i) our results are compatible with those
obtained in [18] (shown in gray): all hybrid EOS fall within this range; ii)
the MFTQCD model allows a deconfinement phase transition at much lower baryon
densities compared to the NJL model (phase transition distribution can be seen
in Fig. 11 in Appendix A). Although, the constraints permitted a phase
transition at a baryon density of the order of the saturation density or
above, the transition does not occur below $\rho_{B}\sim 0.25$ fm-3 within
NJL, while the MFTQCD model admits transitions at saturation density.

Figure 2: Mass-radius diagram of the 90% of CL stiff sets: MFTQCD (left) and
NJL (right). Solid bands represent sets were obtained applying the pQCD
constraint. Observational data is show as: PSR J0030+0451 (cyan and purple)
[56, 57], PSR J0740+6620 (yellow) [58, 59], HESS J1731-347 (magenta) [73], and
GW170817 (light and dark gray) [74].  Figure 3: Corner plots of
$\xi_{\omega\omega}$, the maximum mass ($M_{\text{max}}$) and the speed of
sound at $M_{\text{max}}$ ($c^{2}_{\text{s,max}}$) for the NJL stiff sets. Set
with (without) pQCD constraint is shown in blue (black).

Solving the Tolman-Oppenheimer-Volkov equations [75, 76], we obtain the mass-
radius relation of static spherical and homogeneous neutron stars. Figure 2
shows the mass-radius diagram obtained with the sets built from the stiff
hadronic EOS. The same color codes from Fig. 1 are used. In the mass-radius
diagrams, we also show the 68% confidence zone of the 2D posterior
distribution of PSR J0030+0451 (cyan and magenta) [56, 57] and PSR J0740+6620
(yellow) [58, 59] NICER data, and the HESS J1731-347 (magenta) [73]. The
GW170817 data from the LIGO-Virgo collaboration is represented by the solid
(dashed) gray area with 90% (50%) CL [74].

The pQCD constraint reduces the maximum mass in both NJL and MFTQCD models.
However, this reduction occurs with more intensity in the stiff NJL model.
This effect can be explained by the corner plot given in Fig. 3. In this
figure, we show the relation between $\xi_{\omega\omega}$, the maximum mass
($M_{\text{max}}$) and the speed of sound at the maximum mass
($c^{2}_{s,\text{max}}$) for the NJL stiff set with (without) the pQCD
constraint in blue (black). We identify a strong correlation between these
terms. Thus, $\xi_{\omega\omega}$ is the variable responsible for generating
high values of maximum mass at the cost of considerably increasing the values
of the speed of sound. However, the pQCD constraint is related to causality
[60], allowing only smaller $\xi_{\omega\omega}$. We do not have this behavior
in the soft sets, as the $\xi_{\omega\omega}$ values are smaller (see Figs. 9
and 10 in Appendix A for a corner plot with the model parameters).

The pQCD constraint has also a direct effect on the speed of sound. When this
constraint is imposed, the NJL $c_{s}^{2}$ PDF shows a maximum below
$c_{s}^{2}=0.3$, and the medium value in the centre of the maximum mass
configuration is 0.35, although it can reach values as large as $\sim 0.7$. On
the other hand, for the soft hadron EOS a medium value $\sim 0.7$ is obtained,
with the possibility that a value as large as $\sim 0.85$ is reached in the
centre of the most massive stars (see Table 4). In Table 4, we give several NS
properties (radius, transition density to deconfined quark matter, central
pressure, energy density, baryon density and speed of sound squared) for stars
with a mass equal to 1.4, 1.6, 1.8, 2.0 $M_{\odot}$ and maximum mass.

None of the terms of MFTQCD model affect the speed of sound in a similar way
as the NJL term $\xi_{\omega\omega}$. It is the non-linearity of this term in
the NJL model involving the baryon density that causes the stiffening of the
EOS. MFTQCD EOS predict for $c_{s}^{2}$ in the centre of the most massive
stars values of the order of 0.5, although they can decrease to $\sim 0.35$ in
the case of the stiff EOS.

Figure 4: Comparison between models with pQCD constraint. It is shown, from
the left to the right: the speed of sound as a function of the baryon density,
the mass-radius diagram and the mass by the tidal deformability. The middle
panel with the mass-radius curves also includes observational data as
specified in Fig. 2.

The information on the neutron star properties is completed with Fig. 4 which
shows some properties of all the four sets including the pQCD constraint, soft
and stiff hadron EOS with the NJL and MFTQCD quark models: the speed of sound
squared (left), the mass-radius relations (middle) and the tidal deformability
as a function of the mass (right). The NJL and the MFTQCD sets are shown in
blue and orange, respectively. Sets with the soft (stiff) equations are shown
in solid (patterned) bands. Observational data has the same code color
introduced in Fig. 2. In the tidal deformability versus mass graph (right
panel), the GW170817 data from [74] is represented by the purple line.

Figure 5: 90% CL of the $u,\,d,\,$ and $s$-quark fraction, for MFTQCD (top)
and NJL (bottom) and soft (right) and stiff EOS (right). The bands at the low
density limit define the region where the deconfinement occurs for the
different set EOS, and result from the jump from zero to a finite value of the
order of 1/3 for $u$-quarks and $\lesssim 2/3$ for the $d$-quark.

Some comments are in order: from the Fig. 2 and the Table 4, it is seen that
both NJL and MFTQCD predict similar star radii for masses of the order of
$2\,M_{\odot}$. However, radii below 12 km for medium mass stars can only be
obtained with the MFTQCD model. While within the soft EOS, both quark models
predict for the radius of the maximum mass configuration a value just above 11
km, within the stiff hadron EOS, the NJL predicts almost 13 km, more than 1 km
larger than the MFTQCD model; ii) The pQCD constraints imposes a maximum mass
$\lesssim 2.3\,M_{\odot}$ for MFTQCD and $\lesssim 2.1\,M_{\odot}$ for NJL.
These constraints may lower the maximum mass on about 0.2$M_{\odot}$ for NJL
and 0.06$M_{\odot}$ for MFTQCD for the sets built with the stiff hadronic EOS
and has almost no effect for the soft EOS. See Table 4 where some NS
properties are given for all the sets built in the present work; iii) All sets
are compatible with the NICER observational data. However, only MFTQCD can
partially describe the HESS observation [73]; iv) The NJL stiff set is not
able to describe the tidal deformability obtained from the GW170817 data,
mainly because the stiff hadronic EOS does not satisfy this constraint.
Moreover, the MFTQCD stiff is only partially compatible with it; v) for the
NJL EOS, the deconfiment phase transition occurs just above 2 $\rho_{0}$ and
does not vary much within the complete set. Within the MFTQCD, the
deconfinement phase transition occurs mostly at a lower density, above the
saturation density but below twice the saturation density; vi) for the NJL
EOS, the speed of sound squared suffers a strong decrease at deconfinement,
followed by a steady increase, reaching a value above 0.9 under some
conditions. For the MFTQCD EOS, the increase in the speed of sound that
follows the decrease at the phase transition is weak.

In Fig. 5, the fraction of the three quark flavors is represented as a
function of the density. For NJL, the quark deconfinement occurs before the
onset of the $s$-quark at $\sim 0.4$ (0.5) fm-3 for the stiff (soft) EOS. The
deep in the speed of sound panel for the NJL model after deconfinement results
from the $s$-quark onset. The large baryon density at which is occurs is due
to the fact that the chiral restoration of $s$-quark mass occurs at large
densities. In the MFTQCD model, quarks are already in the chiral restored
phase and it is the larger current mass of the $s$-quark the cause of a
smaller $s$-quark fraction.

Figure 6: The top graphs show the $d_{c}$, and the graphs below show the speed
of sound by baryonic density. In first column, we have model with the soft
equation and in second column, with the stiff equation.

In the following we discuss some more properties of the EOS sets developed in
the present study. In [18, 21, 22], model independent equations were used to
probe quark matter inside NSs. The authors discussed the high density behavior
of the EOS based on the polytropic index $\gamma=\frac{d\ln
P}{d\ln\varepsilon}$, on the normalized matter trace anomaly
$\Delta=\frac{1}{3}-\frac{P}{\varepsilon}$ and some quantities derived from
the trace anomaly. In [22], the authors have conjectured that the matter part
of the trace anomaly could be positive definite and discuss the implication on
the mass-radius curves of NS. In [22], the authors propose that
$d_{c}=\sqrt{\Delta^{2}+\Delta^{\prime 2}}<0.2$ (where $\Delta^{\prime}$ is
$\Delta$’s logarithmic derivative with respect to the energy density
$\varepsilon$) could indicate a deconfined matter behavior. In this work, by
using the microscopic models that define the different datasets, and since the
composition of matter is know, we discuss the implications of these criteria.
Figure 6 shows the 90% CL of $\Delta$ (top panels), $d_{c}$ (middle panels)
and $\gamma$ (bottom panels) as functions of the baryon density represented by
blue (NJL) and orange (MFTQCD) bands. Concerning the trace anomaly, as
discussed in [22], the pQCD constraints favors a positive $\Delta$. Here, we
identify some EOS for which $\Delta$ is always positive. In Fig. 6, the 90% CL
of the $\Delta>0$ at the maximum mass density ($\Delta_{\text{max}}$) EOSs are
represented by dark blue (NJL) and dark orange (MFTQCD) patterns. NJL soft
(stiff) EOS with $\Delta_{\text{max}}>0$ represent 40% (94%) of the respective
sets. For MFTQCD soft (stiff) equations, they represent only 5% (30%). The
high percentage of EOS satisfying $\Delta_{\text{max}}>0$ in the NJL stiff set
is due to the fact that we have low $\rho_{\text{max}}$ values (0.581 to 0.931
fm-3 with 90% CL - see Table 4). In Fig. 6, we see that the band of EOS
without the $\Delta_{\text{max}}>0$ filter already mostly satisfies this
restriction. It is also interesting to note that, when applying the
$\Delta_{\text{max}}>0$ filter, we have a notable change in the phase
transition values of the 90% band in the MFTQCD stiff set. This is because the
MFTQCD stiff set is able to reach higher values of $B$. From Eq. 17, we see
that this parameter directly increases the $\Delta$ values. However, $B$ also
increases the phase transition values, thus creating a correlation between
$\Delta$ and $\rho_{\text{trans}}$. These correlations can be seen in Fig. 8,
which shows the corner plots of the MFTQCD soft (left) and stiff (right) sets.
The effect of $\Delta_{\text{max}}>0$ filter in the mass-radius diagram is
shown in Fig. 7. Patterned (solid) bands show the equations with (without) the
$\Delta_{\text{max}}>0$ filter. In article [21], the authors analyzed the
trace anomaly using agnostic equations. When applying $\Delta>0$, they
observed a decrease in maximum mass and its radius. However, here this
behavior can only be seen in the MFTQCD stiff and NJL soft sets.

Figure 7: Mass-radius diagram of the sets with $\Delta(\rho_{\text{max}})>0$ filter (patterned bands). Sets without this filter are also shown in solid bands. Graphs on the left show the equations with the MFTQCD model, and on the right, with the NJL. |
---|---
Figure 8: Corners of MFTQCD soft (left) and stiff (right). Orange (black)
represents the sets without (with) the $\Delta(\rho_{\text{max}})>0$ filter.

Concerning $d_{c}$ and $\gamma$, we verify that for all sets of models, the
phase transition occurs before the graphs reach the $d_{c}<0.2$ and the
$\gamma<1.75$ region, but rapidly attain a value below these. After crossing
the $d_{c}<0.2$ value, the MFTQCD EOS show a monotonically decreasing behavior
but the same is not true to NJL: the 8-quark term with the
$\xi_{\omega\omega}$ coupling becomes stronger with density and gives rise to
a non-monotonic behavior of $d_{c}$. At large densities $d_{c}$ rises above
0.2. A similar behaviour can be seen in the $\gamma$ values.

NJL
---
| soft w/ pQCD | soft w/o pQCD | stiff w/ pQCD | stiff w/o pQCD
| median | min | max | median | min | max | median | min | max | median | min | max
$M_{\text{max}}$ | 2.024 | 1.893 | 2.108 | 2.023 | 1.893 | 2.109 | 2.058 | 1.968 | 2.130 | 2.120 | 1.979 | 2.311
$\rho_{\text{trans}}$ | 0.356 | 0.315 | 0.388 | 0.356 | 0.315 | 0.389 | 0.375 | 0.324 | 0.400 | 0.371 | 0.315 | 0.398
$R_{\text{max}}$ | 11.26 | 11.09 | 11.75 | 11.26 | 11.09 | 11.75 | 12.97 | 11.79 | 13.71 | 12.31 | 11.61 | 13.64
$P_{0}(M_{\text{max}})$ | 443 | 275 | 539 | 442 | 273 | 540 | 182 | 92 | 390 | 308 | 103 | 513
$\epsilon_{0}(M_{\text{max}})$ | 1272 | 1122 | 1319 | 1271 | 1122 | 1319 | 850 | 640 | 1156 | 1021 | 674 | 1203
$\rho_{0}(M_{\text{max}})$ | 1.015 | 0.942 | 1.041 | 1.014 | 0.942 | 1.041 | 0.733 | 0.581 | 0.931 | 0.833 | 0.607 | 0.945
$c^{2}_{s_{0}}(M_{\text{max}})$ | 0.696 | 0.407 | 0.847 | 0.694 | 0.407 | 0.849 | 0.349 | 0.173 | 0.709 | 0.604 | 0.190 | 0.942
$R(2M_{\odot})$ | 11.75 | 11.35 | 12.01 | 11.75 | 11.34 | 12.00 | 13.25 | 12.39 | 13.75 | 13.43 | 12.67 | 13.73
$P_{0}(2M_{\odot})$ | 221 | 186 | 316 | 220 | 186 | 316 | 99 | 81 | 155 | 97 | 82 | 143
$\epsilon_{0}(2M_{\odot})$ | 851 | 743 | 1081 | 848 | 741 | 1082 | 555 | 390 | 740 | 517 | 396 | 693
$\rho_{0}(2M_{\odot})$ | 0.857 | 0.780 | 0.970 | 0.858 | 0.780 | 0.972 | 0.622 | 0.507 | 0.792 | 0.546 | 0.466 | 0.722
$c^{2}_{s_{0}}(2M_{\odot})$ | 0.609 | 0.544 | 0.732 | 0.608 | 0.544 | 0.735 | 0.335 | 0.223 | 0.539 | 0.370 | 0.253 | 0.482
$R(1.8M_{\odot})$ | 12.67 | 12.48 | 12.76 | 12.67 | 12.48 | 12.77 | 13.76 | 13.50 | 13.76 | 13.76 | 13.47 | 13.76
$P_{0}(1.8M_{\odot})$ | 115 | 108 | 131 | 115 | 108 | 131 | 63 | 63 | 78 | 63 | 63 | 79
$\epsilon_{0}(1.8M_{\odot})$ | 582 | 547 | 670 | 582 | 546 | 669 | 363 | 356 | 470 | 363 | 357 | 465
$\rho_{0}(1.8M_{\odot})$ | 0.548 | 0.519 | 0.618 | 0.547 | 0.518 | 0.617 | 0.359 | 0.353 | 0.450 | 0.359 | 0.354 | 0.446
$c^{2}_{s_{0}}(1.8M_{\odot})$ | 0.336 | 0.266 | 0.383 | 0.335 | 0.266 | 0.383 | 0.404 | 0.209 | 0.501 | 0.363 | 0.232 | 0.498
$R(1.4M_{\odot})$ | 12.98 | 12.89 | 13.00 | 12.97 | 12.89 | 13.00 | 13.57 | 13.57 | 13.57 | 13.57 | 13.57 | 13.57
$P_{0}(1.4M_{\odot})$ | 53 | 52 | 56 | 53 | 52 | 56 | 37 | 37 | 37 | 37 | 37 | 37
$\epsilon_{0}(1.4M_{\odot})$ | 412 | 403 | 426 | 412 | 402 | 426 | 306 | 306 | 306 | 306 | 306 | 308
$\rho_{0}(1.4M_{\odot})$ | 0.407 | 0.399 | 0.420 | 0.407 | 0.398 | 0.420 | 0.309 | 0.309 | 0.309 | 0.309 | 0.309 | 0.310
$c^{2}_{s_{0}}(1.4M_{\odot})$ | 0.341 | 0.288 | 0.361 | 0.341 | 0.288 | 0.361 | 0.396 | 0.275 | 0.409 | 0.392 | 0.260 | 0.408
MFTQCD
| soft w/ pQCD | soft w/o pQCD | stiff w/ pQCD | stiff w/o pQCD
| median | min | max | median | min | max | median | min | max | median | min | max
$M_{\text{max}}$ | 2.101 | 1.968 | 2.218 | 2.105 | 1.957 | 2.212 | 2.143 | 1.955 | 2.318 | 2.161 | 1.965 | 2.378
$\rho_{\text{trans}}$ | 0.259 | 0.173 | 0.365 | 0.256 | 0.174 | 0.369 | 0.316 | 0.185 | 0.388 | 0.311 | 0.183 | 0.388
$R_{\text{max}}$ | 11.27 | 10.78 | 11.61 | 11.27 | 10.77 | 11.60 | 11.80 | 11.05 | 13.32 | 11.86 | 11.07 | 12.74
$P_{0}(M_{\text{max}})$ | 432 | 408 | 459 | 432 | 410 | 459 | 396 | 88 | 445 | 393 | 89 | 444
$\epsilon_{0}(M_{\text{max}})$ | 1239 | 1135 | 1370 | 1237 | 1140 | 1378 | 1135 | 791 | 1312 | 1126 | 802 | 1303
$\rho_{0}(M_{\text{max}})$ | 0.993 | 0.915 | 1.100 | 0.991 | 0.919 | 1.104 | 0.913 | 0.706 | 1.050 | 0.906 | 0.712 | 1.043
$c^{2}_{s_{0}}(M_{\text{max}})$ | 0.488 | 0.465 | 0.503 | 0.488 | 0.463 | 0.503 | 0.495 | 0.355 | 0.521 | 0.498 | 0.366 | 0.528
$R(2M_{\odot})$ | 12.01 | 11.43 | 12.37 | 12.00 | 11.46 | 12.35 | 12.68 | 11.72 | 13.49 | 12.77 | 11.78 | 13.57
$P_{0}(2M_{\odot})$ | 199 | 148 | 339 | 195 | 150 | 321 | 145 | 86 | 275 | 138 | 86 | 272
$\epsilon_{0}(2M_{\odot})$ | 742 | 598 | 1086 | 732 | 602 | 1047 | 621 | 483 | 947 | 605 | 447 | 940
$\rho_{0}(2M_{\odot})$ | 0.714 | 0.601 | 0.908 | 0.711 | 0.609 | 0.897 | 0.613 | 0.502 | 0.844 | 0.586 | 0.464 | 0.830
$c^{2}_{s_{0}}(2M_{\odot})$ | 0.471 | 0.465 | 0.479 | 0.471 | 0.465 | 0.478 | 0.473 | 0.463 | 0.478 | 0.473 | 0.464 | 0.481
$R(1.8M_{\odot})$ | 12.28 | 11.64 | 12.64 | 12.28 | 11.61 | 12.66 | 13.17 | 12.00 | 13.76 | 13.19 | 12.03 | 13.76
$P_{0}(1.8M_{\odot})$ | 126 | 103 | 175 | 125 | 104 | 179 | 88 | 63 | 148 | 86 | 63 | 145
$\epsilon_{0}(1.8M_{\odot})$ | 586 | 499 | 742 | 583 | 503 | 756 | 483 | 362 | 668 | 475 | 363 | 661
$\rho_{0}(1.8M_{\odot})$ | 0.555 | 0.484 | 0.683 | 0.553 | 0.487 | 0.693 | 0.465 | 0.358 | 0.621 | 0.457 | 0.358 | 0.615
$c^{2}_{s_{0}}(1.8M_{\odot})$ | 0.443 | 0.432 | 0.451 | 0.443 | 0.430 | 0.451 | 0.447 | 0.363 | 0.497 | 0.449 | 0.385 | 0.497
$R(1.4M_{\odot})$ | 12.42 | 11.81 | 12.93 | 12.41 | 11.78 | 12.94 | 13.57 | 12.09 | 13.57 | 13.57 | 12.20 | 13.57
$P_{0}(1.4M_{\odot})$ | 64 | 55 | 79 | 63 | 55 | 80 | 37 | 37 | 69 | 38 | 37 | 68
$\epsilon_{0}(1.4M_{\odot})$ | 441 | 399 | 508 | 441 | 400 | 511 | 307 | 306 | 467 | 329 | 306 | 468
$\rho_{0}(1.4M_{\odot})$ | 0.437 | 0.400 | 0.498 | 0.436 | 0.401 | 0.502 | 0.310 | 0.309 | 0.461 | 0.332 | 0.309 | 0.460
$c^{2}_{s_{0}}(1.4M_{\odot})$ | 0.426 | 0.407 | 0.437 | 0.426 | 0.407 | 0.437 | 0.399 | 0.309 | 0.447 | 0.401 | 0.295 | 0.451
Table 4: The median and the 90% CL of NS properties for all sets: maximum mass
$M_{\text{max}}$ ($M_{\odot}$); transition baryon density
$\rho_{\text{trans}}$ (fm-3); radium $R_{\text{i}}$ (km), central pressure
$P_{0,i}$ (MeV/fm-3), central energy density $\epsilon_{0,i}$ (MeV/fm-3),
central baryonic density $\rho_{0,i}$ (fm-3) at $M_{i}$, with
$i=[\text{max},2M_{\odot},1.8M_{\odot},1.4M_{\odot}]$.

##  V Conclusions

In the present study we have as main objective to understand the possibility
that NS have a deconfined quark phase, considering the quark phase described
by phenomenological microscopic models, and avoiding agnostic descriptions.
For the hadron phase, we use two extreme hadron EOS from [25]: the BMPF 220
(soft) and the BMPF 260 (stiff) that satisfy nuclear matter properties and
observational constraints. For the quark phase, we consider the chiral
symmetric SU(3) Nambu-Jona-Lasinio (NJL) model with 4-quark and 8-quark
interaction terms, and the Mean Field Theory of QCD (MFTQCD) which results
from the QCD Lagrangian, considering a gluon field decomposition in soft and
hard momentum components [46, 63, 64].

Based on these models we have generated eight different sets choosing two
hadron EOS and the two quark model, and considering/not considering the very
high density pQCD constraints. Bayesian inference is applied to determine the
quark models coupling constants and a constant bag term $B$. We have imposed
as constraints the observational data from PSR J0030+0451 [56, 57] and PSR
J0740+6620 [58, 59]. Besides, we have also imposed that the phase transition
should occur around $\rho_{\text{trans}}\approx 0.15-0.4\text{ fm}^{-3}$. In
four sets, we have also enforced the pQCD constraint at $\rho_{B}=7n_{0}$
(where $n_{0}=0.16$ fm-3) [60]. We have confirmed that the current
observational data are compatible with the existence of a quark core inside
NS. The pQCD constraints were important to constrain the NJL model to give
causal EOS, since the NJL 8-quark term which is responsible for allowing two
solar mass stars with a large quark core may give rise to a super-luminous
speed of sound. For both quark models we have obtained for the maximum star
mass a value of the order of 2.1–2.3$M_{\odot}$ and similar radii. However,
with MFTQCD quark model it was possible to obtain medium and low mass NS with
a smaller radius, in particular, below 12 km (even compatible with the compact
object HESS J1731-347 [73]). With MTFQCD we have obtained a transition to
quark matter above $\sim 0.15$ fm-3 but within NJL the transition density was
above $\sim 0.25$ fm-3. This explains why within NJL the medium and low mass
stars have a larger radius. Note that in [38], only four quark interactions
were introduced in the NJL Lagrangian density, and, as a consequence, quark
cores were only found in the most massive stars. In [31], multiquark
interactions were considered in the NJL model, however, the main objective of
the study was to determine under which conditions a NS ”third family” with a
disconnected hybrid star branch exists. In our analysis, the multiquark
interactions considered do not predict a third family, since we did not
include a eight quark scalar-pseudoscalar interaction because we did not want
to affect the vacuum properties.

We have discussed the behavior of the normalized matter trace anomaly and the
derived quantity $d_{c}$ proposed to measure conformality in [22] and could
conclude that NS may have a quark core but within the present description,
quark matter is still strongly interacting matter and the conformal limit is
not attained in the NS center. Besides the multiquark vector interactions may
give rise to values of $d_{c}>0.2$ and $\gamma>1.7$ in a deconfined phase.

In a future work, we will apply the Bayesian inference also in the hadronic
phase. In that way, we will obtain all the range of the equations of state. We
also plan to compare to purely hadronic sets and analyze if hybrid stars are
more probable than hadronic stars based on the current observational data.

## Acknowledgments

M.A. expresses sincere gratitude to the FCT for their generous support through
Ph.D. grant number 2022.11685.BD. This research received partial funding from
national sources through FCT (Fundação para a Ciência e a Tecnologia, I.P,
Portugal) for projects UIDB/04564/2020 and UIDP/04564/2020, identified by DOI
10.54499/UIDB/04564/2020 and 10.54499/UIDP/04564/2020, respectively, as well
as for project 2022.06460.PTDC with DOI 10.54499/2022.06460.PTDC. The authors
acknowledge the Laboratory for Advanced Computing at the University of Coimbra
for providing HPC resources that have contributed to the research results
reported within this paper, URL: https://www.uc.pt/lca.

##  Appendix A Parameter and phase transition posteriors

In this section, we present some additional results that support our
discussion. Figures 9 and 10 show the corner plots of the NJL and MFTQCD
models parameters, respectively. Each corner shows a set with the pQCD
constraint (blue or orange) and the same set without the constraint (black).
Figure 11 shows the probability distribution of the baryonic density of phase
transition. Plots are arranged as shown in figure 9 and 10.

|
---|---
Figure 9: Corner plots of the parameters values of NJL soft (left) and NJL stiff (right). Models with pQCD constraints are plotted in blue and without pQCD in black. |
---|---
Figure 10: Corner plots of the parameters values of MFTQCD soft (left) and MFTQCD stiff (right). Models with pQCD constraints are plotted in orange and without pQCD in black. |
---|---
|
Figure 11: The phase transition likelihood ($\mathcal{L}^{\rm PhT}$) and the
marginalized posteriors. The upper figures show the graphs for the NJL model
and the lower ones for the MFTQCD model. On the left, we have the soft models
and on the right, the stiff models. Models with pQCD constraints are plotted
in blue (NJL) or orange (MFTQCD) and without pQCD in black. Dotted lines shows
the imposed constraint in phase transition ($\mathcal{L}^{\text{PhT}}$).

## References

  * Fukushima and Hatsuda [2011] K. Fukushima and T. Hatsuda, The phase diagram of dense QCD, Rept. Prog. Phys. 74, 014001 (2011), arXiv:1005.4814 [hep-ph] .
  * Bazavov  _et al._ [2012] A. Bazavov _et al._ , The chiral and deconfinement aspects of the QCD transition, Phys. Rev. D 85, 054503 (2012), arXiv:1111.1710 [hep-lat] .
  * Glendenning [1997] N. K. Glendenning,  _Compact stars: Nuclear physics, particle physics, and general relativity_ (1997).
  * Rezzolla  _et al._ [2018] L. Rezzolla, P. Pizzochero, D. I. Jones, N. Rea, and I. Vidaña, eds.,  _The Physics and Astrophysics of Neutron Stars_, Vol. 457 (Springer, 2018).
  * Kurkela  _et al._ [2010] A. Kurkela, P. Romatschke, and A. Vuorinen, Cold Quark Matter, Phys. Rev. D 81, 105021 (2010), arXiv:0912.1856 [hep-ph] .
  * Hebeler  _et al._ [2013] K. Hebeler, J. M. Lattimer, C. J. Pethick, and A. Schwenk, Equation of state and neutron star properties constrained by nuclear physics and observation, Astrophys. J. 773, 11 (2013), arXiv:1303.4662 [astro-ph.SR] .
  * Drischler  _et al._ [2016] C. Drischler, K. Hebeler, and A. Schwenk, Asymmetric nuclear matter based on chiral two- and three-nucleon interactions, Phys. Rev. C 93, 054314 (2016), arXiv:1510.06728 [nucl-th] .
  * Son and Stephanov [2001] D. T. Son and M. A. Stephanov, QCD at finite isospin density, Phys. Rev. Lett. 86, 592 (2001), arXiv:hep-ph/0005225 .
  * Ebert and Klimenko [2006] D. Ebert and K. G. Klimenko, Pion condensation in electrically neutral cold matter with finite baryon density, Eur. Phys. J. C 46, 771 (2006), arXiv:hep-ph/0510222 .
  * Barducci  _et al._ [2004] A. Barducci, R. Casalbuoni, G. Pettini, and L. Ravagli, A Calculation of the QCD phase diagram at finite temperature, and baryon and isospin chemical potentials, Phys. Rev. D 69, 096004 (2004), arXiv:hep-ph/0402104 .
  * Alford  _et al._ [1998] M. G. Alford, K. Rajagopal, and F. Wilczek, QCD at finite baryon density: Nucleon droplets and color superconductivity, Phys. Lett. B 422, 247 (1998), arXiv:hep-ph/9711395 .
  * Mishra and Mishra [2004] A. Mishra and H. Mishra, Chiral symmetry breaking, color superconductivity and color neutral quark matter: A Variational approach, Phys. Rev. D 69, 014014 (2004), arXiv:hep-ph/0306105 .
  * Bonanno and Sedrakian [2012] L. Bonanno and A. Sedrakian, Composition and stability of hybrid stars with hyperons and quark color-superconductivity, Astron. Astrophys. 539, A16 (2012), arXiv:1108.0559 [astro-ph.SR] .
  * Abhishek and Mishra [2021] A. Abhishek and H. Mishra, Chiral Symmetry Breaking, Color Superconductivity, and Equation of State for Magnetized Strange Quark Matter, Springer Proc. Phys. 261, 593 (2021).
  * Alford  _et al._ [1999] M. G. Alford, K. Rajagopal, and F. Wilczek, Color flavor locking and chiral symmetry breaking in high density QCD, Nucl. Phys. B 537, 443 (1999), arXiv:hep-ph/9804403 .
  * Mannarelli  _et al._ [2006] M. Mannarelli, K. Rajagopal, and R. Sharma, Testing the Ginzburg-Landau approximation for three-flavor crystalline color superconductivity, Phys. Rev. D 73, 114012 (2006), arXiv:hep-ph/0603076 .
  * Rajagopal and Sharma [2006] K. Rajagopal and R. Sharma, The Crystallography of Three-Flavor Quark Matter, Phys. Rev. D 74, 094019 (2006), arXiv:hep-ph/0605316 .
  * Annala  _et al._ [2020] E. Annala, T. Gorda, A. Kurkela, J. Nättilä, and A. Vuorinen, Evidence for quark-matter cores in massive neutron stars, Nature Phys. 16, 907 (2020), arXiv:1903.09121 [astro-ph.HE] .
  * Altiparmak  _et al._ [2022] S. Altiparmak, C. Ecker, and L. Rezzolla, On the Sound Speed in Neutron Stars, Astrophys. J. Lett. 939, L34 (2022), arXiv:2203.14974 [astro-ph.HE] .
  * Somasundaram  _et al._ [2023] R. Somasundaram, I. Tews, and J. Margueron, Perturbative QCD and the neutron star equation of state, Phys. Rev. C 107, L052801 (2023), arXiv:2204.14039 [nucl-th] .
  * Fujimoto  _et al._ [2022] Y. Fujimoto, K. Fukushima, L. D. McLerran, and M. Praszalowicz, Trace Anomaly as Signature of Conformality in Neutron Stars, Phys. Rev. Lett. 129, 252702 (2022), arXiv:2207.06753 [nucl-th] .
  * Annala  _et al._ [2023] E. Annala, T. Gorda, J. Hirvonen, O. Komoltsev, A. Kurkela, J. Nättilä, and A. Vuorinen, Strongly interacting matter exhibits deconfined behavior in massive neutron stars, Nature Commun. 14, 8451 (2023), arXiv:2303.11356 [astro-ph.HE] .
  * Baym and Chin [1976] G. Baym and S. Chin, Can a neutron star be a giant mit bag?, Physics Letters B 62, 241 (1976).
  * Oertel  _et al._ [2017] M. Oertel, M. Hempel, T. Klähn, and S. Typel, Equations of state for supernovae and compact stars, Rev. Mod. Phys. 89, 015007 (2017), arXiv:1610.03361 [astro-ph.HE] .
  * Malik  _et al._ [2023] T. Malik, M. Ferreira, M. B. Albino, and C. Providência, Spanning the full range of neutron star properties within a microscopic description, Physical Review D 107, 10.1103/physrevd.107.103018 (2023).
  * Malik  _et al._ [2022] T. Malik, M. Ferreira, B. K. Agrawal, and C. Providência, Relativistic description of dense matter equation of state and compatibility with neutron star observables: A bayesian approach, The Astrophysical Journal 930, 17 (2022).
  * Logoteta  _et al._ [2013] D. Logoteta, C. Providência, and I. Vidaña, Formation of hybrid stars from metastable hadronic stars, Phys. Rev. C 88, 055802 (2013), arXiv:1311.0618 [nucl-th] .
  * Benic  _et al._ [2015] S. Benic, D. Blaschke, D. E. Alvarez-Castillo, T. Fischer, and S. Typel, A new quark-hadron hybrid equation of state for astrophysics - I. High-mass twin compact stars, Astron. Astrophys. 577, A40 (2015), arXiv:1411.2856 [astro-ph.HE] .
  * Ferreira  _et al._ [2021] M. Ferreira, R. C. Pereira, and C. Providência, Hybrid stars with large strange quark cores constrained by GW170817, Physical Review D 103, 10.1103/physrevd.103.123020 (2021).
  * Matsuoka  _et al._ [2018] H. Matsuoka, Y. Tsue, J. a. Da Providência, C. Providência, and M. Yamamura, Hybrid stars from the NJL model with a tensor-interaction, Phys. Rev. D 98, 074027 (2018), arXiv:1806.04377 [hep-ph] .
  * Alvarez-Castillo  _et al._ [2016] D. Alvarez-Castillo, A. Ayriyan, S. Benic, D. Blaschke, H. Grigorian, and S. Typel, New class of hybrid EoS and Bayesian M-R data analysis, Eur. Phys. J. A 52, 69 (2016), arXiv:1603.03457 [nucl-th] .
  * Han and Steiner [2019] S. Han and A. W. Steiner, Tidal deformability with sharp phase transitions in binary neutron stars, Phys. Rev. D 99, 083014 (2019).
  * Schramm  _et al._ [2016] S. Schramm, V. Dexheimer, and R. Negreiros, Modelling hybrid stars in quark-hadron approaches, The European Physical Journal A 52, 10.1140/epja/i2016-16014-5 (2016).
  * Parisi  _et al._ [2020] A. Parisi, C. V. Flores, C. H. Lenzi, C.-S. Chen, and G. Lugones, Hybrid stars in the light of the merging event gw170817 (2020), arXiv:2009.14274 [astro-ph.HE] .
  * Tonetto and Lugones [2020] L. Tonetto and G. Lugones, Discontinuity gravity modes in hybrid stars: Assessing the role of rapid and slow phase conversions, Phys. Rev. D 101, 123029 (2020).
  * Bastian [2021] N.-U. F. Bastian, Phenomenological quark-hadron equations of state with first-order phase transitions for astrophysical applications, Phys. Rev. D 103, 023001 (2021).
  * Minamikawa  _et al._ [2021] T. Minamikawa, T. Kojo, and M. Harada, Quark-hadron crossover equations of state for neutron stars: Constraining the chiral invariant mass in a parity doublet model, Phys. Rev. C 103, 045205 (2021).
  * Pfaff  _et al._ [2022] A. Pfaff, H. Hansen, and F. Gulminelli, Bayesian analysis of the properties of hybrid stars with the Nambu–Jona-Lasinio model, Phys. Rev. C 105, 035802 (2022), arXiv:2112.09595 [astro-ph.HE] .
  * Xie and Li [2021] W.-J. Xie and B.-A. Li, Bayesian inference of the dense-matter equation of state encapsulating a first-order hadron-quark phase transition from observables of canonical neutron stars, Phys. Rev. C 103, 035802 (2021), arXiv:2009.13653 [nucl-th] .
  * Ayriyan  _et al._ [2021] A. Ayriyan, D. Blaschke, A. G. Grunfeld, D. Alvarez-Castillo, H. Grigorian, and V. Abgaryan, Bayesian analysis of multimessenger M-R data with interpolated hybrid EoS, Eur. Phys. J. A 57, 318 (2021), arXiv:2102.13485 [astro-ph.HE] .
  * Takatsy  _et al._ [2023] J. Takatsy, P. Kovacs, G. Wolf, and J. Schaffner-Bielich, What neutron stars tell about the hadron-quark phase transition: A Bayesian study, Phys. Rev. D 108, 043002 (2023), arXiv:2303.00013 [astro-ph.HE] .
  * Klevansky [1992] S. P. Klevansky, The Nambu-Jona-Lasinio model of quantum chromodynamics, Rev. Mod. Phys. 64, 649 (1992).
  * Hatsuda and Kunihiro [1994] T. Hatsuda and T. Kunihiro, QCD phenomenology based on a chiral effective Lagrangian, Phys. Rept. 247, 221 (1994), arXiv:hep-ph/9401310 .
  * Buballa and Oertel [1999] M. Buballa and M. Oertel, Strange quark matter with dynamically generated quark masses, Phys. Lett. B 457, 261 (1999), arXiv:hep-ph/9810529 .
  * Ferreira  _et al._ [2020a] M. Ferreira, R. Câmara Pereira, and C. Providência, Quark matter in light neutron stars, Phys. Rev. D 102, 083030 (2020a), arXiv:2008.12563 [nucl-th] .
  * Fogaca and Navarra [2011] D. A. Fogaca and F. S. Navarra, Gluon condensates in a cold quark–gluon plasma, Phys. Lett. B 700, 236 (2011), arXiv:1012.5266 [hep-ph] .
  * Nambu and Jona-Lasinio [1961a] Y. Nambu and G. Jona-Lasinio, Dynamical Model of Elementary Particles Based on an Analogy with Superconductivity. 1., Phys. Rev. 122, 345 (1961a).
  * Nambu and Jona-Lasinio [1961b] Y. Nambu and G. Jona-Lasinio, Dynamical model of elementary particles based on an analogy with superconductivity. II., Phys. Rev. 124, 246 (1961b).
  * Buballa [2005a] M. Buballa, NJL model analysis of quark matter at large density, Phys. Rept. 407, 205 (2005a), arXiv:hep-ph/0402234 .
  * Schertler  _et al._ [1999] K. Schertler, S. Leupold, and J. Schaffner-Bielich, Neutron stars and quark phases in the NJL model, Phys. Rev. C 60, 025801 (1999), arXiv:astro-ph/9901152 .
  * Baldo  _et al._ [2003] M. Baldo, M. Buballa, F. Burgio, F. Neumann, M. Oertel, and H. J. Schulze, Neutron stars and the transition to color superconducting quark matter, Phys. Lett. B 562, 153 (2003), arXiv:nucl-th/0212096 .
  * Menezes and Providencia [2003] D. P. Menezes and C. Providencia, Warm stellar matter with deconfinement: Application to compact stars, Phys. Rev. C 68, 035804 (2003), arXiv:nucl-th/0308041 .
  * Shovkovy  _et al._ [2003] I. Shovkovy, M. Hanauske, and M. Huang, Nonstrange hybrid compact stars with color superconducting matter, Phys. Rev. D 67, 103004 (2003), arXiv:hep-ph/0303027 .
  * Pagliara and Schaffner-Bielich [2008] G. Pagliara and J. Schaffner-Bielich, Stability of CFL cores in Hybrid Stars, Phys. Rev. D 77, 063004 (2008), arXiv:0711.1119 [astro-ph] .
  * Ferreira  _et al._ [2020b] M. Ferreira, R. Câmara Pereira, and C. Providência, Neutron stars with large quark cores, Phys. Rev. D 101, 123030 (2020b), arXiv:2005.10543 [nucl-th] .
  * Riley  _et al._ [2019] T. E. Riley _et al._ , A $NICER$ View of PSR J0030+0451: Millisecond Pulsar Parameter Estimation, Astrophys. J. Lett. 887, L21 (2019), arXiv:1912.05702 [astro-ph.HE] .
  * Miller  _et al._ [2019] M. C. Miller _et al._ , PSR J0030+0451 Mass and Radius from $NICER$ Data and Implications for the Properties of Neutron Star Matter, Astrophys. J. Lett. 887, L24 (2019), arXiv:1912.05705 [astro-ph.HE] .
  * Riley  _et al._ [2021] T. E. Riley _et al._ , A NICER View of the Massive Pulsar PSR J0740+6620 Informed by Radio Timing and XMM-Newton Spectroscopy, Astrophys. J. Lett. 918, L27 (2021), arXiv:2105.06980 [astro-ph.HE] .
  * Miller  _et al._ [2021] M. C. Miller _et al._ , The Radius of PSR J0740+6620 from NICER and XMM-Newton Data, Astrophys. J. Lett. 918, L28 (2021), arXiv:2105.06979 [astro-ph.HE] .
  * Komoltsev and Kurkela [2022] O. Komoltsev and A. Kurkela, How Perturbative QCD Constrains the Equation of State at Neutron-Star Densities, Phys. Rev. Lett. 128, 202701 (2022), arXiv:2111.05350 [nucl-th] .
  * Buballa [2005b] M. Buballa, NJL-model analysis of dense quark matter, Physics Reports 407, 205 (2005b).
  * Olive [2014] K. Olive, Review of particle physics, Chinese Physics C 38, 090001 (2014).
  * Franzon  _et al._ [2012] B. Franzon, D. A. Fogaca, F. S. Navarra, and J. E. Horvath, Self-bound Interacting QCD Matter in Compact Stars, Phys. Rev. D 86, 065031 (2012), arXiv:1203.6090 [astro-ph.SR] .
  * Albino  _et al._ [2021] M. B. Albino, R. Fariello, and F. S. Navarra, Tidal Deformability of Quark Stars with Repulsive Interactions, Phys. Rev. D 104, 083011 (2021), arXiv:2106.12956 [nucl-th] .
  * Celenza and Shakin [1986] L. S. Celenza and C. M. Shakin, Description of the gluon condensate, Phys. Rev. D 34, 1591 (1986).
  * Li and Shakin [2005] X. Li and C. M. Shakin, Description of gluon propagation in the presence of an ${A}^{2}$ condensate, Phys. Rev. D 71, 074007 (2005).
  * Serot and Walecka [1986] B. D. Serot and J. D. Walecka, The Relativistic Nuclear Many Body Problem, Adv. Nucl. Phys. 16, 1 (1986).
  * Tezuka [1987] H. Tezuka, Mean Field Approximation to QCD. 1. Selfconsistent equations of motion, (1987).
  * Landry  _et al._ [2020] P. Landry, R. Essick, and K. Chatziioannou, Nonparametric constraints on neutron star matter with existing and upcoming gravitational wave and pulsar observations, Phys. Rev. D 101, 123007 (2020), arXiv:2003.04880 [astro-ph.HE] .
  * Buchner  _et al._ [2014] J. Buchner, A. Georgakakis, K. Nandra, L. Hsu, C. Rangel, M. Brightman, A. Merloni, M. Salvato, J. Donley, and D. Kocevski, X-ray spectral modelling of the AGN obscuring region in the CDFS: Bayesian model selection and catalogue, Astronomy & Astrophysics 564, A125 (2014), arXiv:1402.0004 [astro-ph.HE] .
  * Buchner [2023] J. Buchner, Nested sampling methods, Statistics Surveys 17, 10.1214/23-ss144 (2023).
  * Ashton  _et al._ [2019] G. Ashton, M. Hübner, P. D. Lasky, C. Talbot, K. Ackley, S. Biscoveanu, Q. Chu, A. Divakarla, P. J. Easter, B. Goncharov, F. H. Vivanco, J. Harms, M. E. Lower, G. D. Meadors, D. Melchor, E. Payne, M. D. Pitkin, J. Powell, N. Sarin, R. J. E. Smith, and E. Thrane, Bilby: A user-friendly bayesian inference library for gravitational-wave astronomy, The Astrophysical Journal Supplement Series 241, 27 (2019).
  * Doroshenko  _et al._ [2022] V. Doroshenko, V. Suleimanov, G. Pühlhofer, and A. Santangelo, A strangely light neutron star within a supernova remnant, Nature Astronomy 6, 1444 (2022).
  * Abbott  _et al._ [2018] B. P. Abbott _et al._ (LIGO Scientific, Virgo), GW170817: Measurements of neutron star radii and equation of state, Phys. Rev. Lett. 121, 161101 (2018), arXiv:1805.11581 [gr-qc] .
  * Tolman [1939] R. C. Tolman, Static solutions of Einstein’s field equations for spheres of fluid, Phys. Rev. 55, 364 (1939).
  * Oppenheimer and Volkoff [1939] J. R. Oppenheimer and G. M. Volkoff, On massive neutron cores, Phys. Rev. 55, 374 (1939).