A Bayesian Inference of Hybrid Stars with Large Quark Cores

Abstract

Neutron stars (NSs) are interesting objects capable of reaching densities unattainable on Earth. The properties of matter under these conditions remain a mystery. Exotic matter, including quark matter, may be present in the NS core. In this work, we explore the possible compositions of NS cores, in particular, the possible existence of large quark cores. We use the Relativistic Mean Field (RMF) model with nonlinear terms for the hadron phase and the Nambu–Jona-Lasinio (NJL) model and Mean Field Theory of Quantum Chromodynamics (MFTQCD) for the quark phase. Through Bayesian inference, we obtain different sets of equations: four sets with hybrid equations (three using the NJL model and the other using the MFTQCD model), and one set with only the hadron phase. We impose constraints regarding the properties of nuclear matter, X-ray observational data from NICER, perturbative QCD (pQCD) calculations, and causality on all sets. One set of hybrid NJL equations of state was also constrained by adding the GW170817 detection. All sets can describe observational data and theoretical restrictions. The MFTQCD allows for a phase transition to quark matter at lower densities compared to the NJL models. The MFTQCD model indicates that NSs with 1.4 $\text{M}_{\odot}$ have quark matter in their inner core. However, NJL models suggest that it is more probable that 1.4 $\text{M}_{\odot}$ NSs do not contain quark matter. Both the MFTQCD and NJL models agree that there is quark matter in 2 $\text{M}_{\odot}$ NSs. It is discussed that hybrid stars with a stiff quark equation of state could explain a larger radius of more massive stars, such as two solar mass stars, with respect to the canonical NS.

I Introduction

The internal composition of neutron stars (NSs) remains one of the most significant open questions in nuclear astrophysics. These compact objects, with masses of approximately 1.2–2.0 solar masses ($\text{M}_{\odot}$) concentrated within a radius of only 10-14 km, represent the densest observable matter in the universe [1, 2]. At such extreme densities, exceeding several times the nuclear saturation density ($\rho_{\text{sat}}\approx 2.7\times 10^{14}$ g cm^-3) our understanding of matter’s behavior becomes increasingly uncertain due to the limitations of terrestrial experiments and first-principles calculations. NSs are believed to provide a natural laboratory for studying exotic high-density phases of QCD, such as the neutron superfluid phase [2, 3]. The outer cores of NSs are so dense and thick that electromagnetic signals cannot escape, while theoretical calculations of the inner cores are hindered by the limitations of first principles lattice Quantum Chromodynamics (QCD).

A particularly intriguing possibility is that NSs may undergo a phase transition from hadronic matter to quark matter in their inner cores, forming what are known as hybrid stars [1, 4, 5, 6, 7, 8]. Quantum Chromodynamics (QCD), the fundamental theory of strong interactions, predicts a transition from confined hadronic matter to a deconfined quark-gluon plasma at sufficiently high densities or temperatures. While numerical simulations of QCD at vanishing baryonic chemical potential indicate a smooth crossover transition at a temperature of $T\approx 154.9$ MeV [9], the nature of this transition at the high densities and relatively low temperatures relevant for NS interiors remains an open question. Some studies suggest that finite surface tension effects can lead to mixed phase states with different geometric shapes (known as “pasta” phases), potentially inducing a smooth phase transition [10, 11, 12]. Depending on the nature of the phase transition, a third family of stable, compact stars (twin stars) with different radii compared to normal NSs may appear, providing a unique observational signature of the hadron-quark transition [13, 12].

Recent breakthroughs in multi-messenger astronomy have opened unprecedented opportunities to probe the properties of supranuclear matter. The detection of gravitational waves (GWs) from binary NS mergers, beginning with GW170817 [14, 15, 16], coupled with electromagnetic observations [17], has provided valuable constraints on the NS equation of state (EOS). Additionally, precise mass and radius measurements from NASA’s Neutron Star Interior Composition Explorer (NICER) mission have further constrained the possible EOS models [18, 19, 12]. The accuracy of pulsar timing, comparable to that of atomic clocks (one part in 10¹⁵), allows for the indirect detection of GWs from binary NS merger events and provides a window for exploring phase transitions occurring inside a pulsar core [20]. Notably, the GW signal emitted during the final orbits of colliding NSs contains imprints of the tidal deformability parameter $\Lambda$, which can be related to the properties of dense matter in terms of the EOS [21, 22, 23].

These observational advances have motivated a renewed theoretical effort to develop more sophisticated EOS models that can account for potential phase transitions while remaining consistent with observational constraints. In this context, Bayesian inference has emerged as a powerful framework for parameter estimation and model selection in astrophysics [24]. It provides a natural way to incorporate prior knowledge, handle uncertainties, and update our beliefs based on new observations.

In this work, we present a comprehensive Bayesian analysis of the EOS for hybrid stars, simultaneously sampling the parameters of both the hadron and quark matter phases. For the hadron phase, we employ the Relativistic Mean Field (RMF) theory with non-linear terms, which has been widely used to describe nuclear matter and finite nuclei [25, 26, 27, 28, 29, 30, 31]. This approach uses the RMF model to describe the hadron phase of NS matter, involving baryons interacting through the exchange of mesons. For the quark matter phase, we adopt two models: the Nambu–Jona-Lasinio (NJL) model [32, 33, 34, 35] and the Mean Field Theory of QCD (MFTQCD) model [36]. The first one incorporates key features of QCD such as global symmetries from QCD and dynamical chiral symmetry breaking and its manifestations [35, 37, 38, 39]. The NJL model offers an attractive framework as it describes the quark phase with a three-flavor version and parameters determined by fitting to various meson and baryon masses [39, 40]. The second model is obtained by decomposing the gluon field into low- and high-momentum components and applying suitable approximations to these fields. The MFTQCD EOS exhibits behavior similar to that of the vectorial MIT bag model [41] and can result in an mass-radius diagram consistent with recent observations [42, 43]. This combined approach allows us to explore a more complete picture of the dense matter in NS cores.

The existence of phase transitions in NS cores can manifest in observable signatures across multiple messengers. In the mass-radius diagram, a strong first-order phase transition can produce disconnected branches of stable configurations, leading to the possibility of “twin stars” - NSs with the same mass but different radii [13, 12]. During binary NS mergers, the hadron-quark phase transition can significantly affect the dynamics of the system and the emitted GWs [44, 45].

The transition from hadron to quark matter in our model employs the Maxwell construction for phase equilibrium, where the transition occurs at a specific pressure with a discontinuity in energy density [1, 46]. While the Gibbs construction allows for a mixed phase where hadronic and quark matter coexist, the Maxwell construction assumes a sharp interface between the two phases with equal pressures but different densities. Our choice of the Maxwell construction is motivated by several considerations. First, it provides a more conservative estimate of the transition effects, as the energy density discontinuity leads to more pronounced observational signatures in GW emission [47]. Second, the surface tension at the hadron-quark interface, though poorly constrained, is believed to be sufficiently large to disfavor the formation of mixed-phase structures in many scenarios [11, 48]. Third, the simplified thermodynamic treatment of the Maxwell construction is computationally advantageous for our Bayesian parameter estimation, allowing for more extensive sampling of the parameter space.

By performing this comprehensive Bayesian inference of the hybrid star EOS, we aim to address several key questions: (1) Is it possible to build hybrid EOSs that satisfy current observations and theoretical predictions while containing a large quark core? (2) What are the most likely properties of this transition, such as its onset density and strength? (3) Which properties distinguish purely hadronic from hybrid EOS models?

The article is organized as follows. In Sec. II, we describe the theoretical framework for the RMF, NJL and MFTQCD models and the construction of the hybrid EOS. Sec. III outlines our Bayesian methodology, including the prior distributions, likelihood function, and computational techniques. In Sec. IV, we present our results on the posterior distributions of model parameters and the resulting constraints on the hybrid EOS. Sec. V discusses the implications of our findings for NS observations and fundamental nuclear physics. Finally, Sec. VI summarizes our conclusions and outlines directions for future work.

II Equations of State

This section describes the models used for the hadron and quark phases. Hybrid EOSs are built through Maxwell construction. Four sets of hybrid EOSs were generated: three using the NJL model for the quark phase with different priors and constraints and one using the MFTQCD for the quark phase. A fifth set consisting only of nucleonic matter was also generated. The hadron phase is described by the RMF model in all cases.

II.1 Hadron phase

The Relativistic Mean Field (RMF) model is used for the hadron phase. In this model, nucleon interactions are mediated by the exchange of mesons: the scalar-isoscalar meson $\sigma$, the vector-isoscalar meson $\omega$, and the vector-isovector meson $\varrho$. This work also adds non-linear terms. The Lagrangian is given by [31]:

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

where

$$\mathcal{L}_{N} =\bar{\Psi}\left[\gamma^{\mu}\left(i\partial_{\mu}-g_{\omega}\omega_{\mu}-g_{\rho}\boldsymbol{t}\cdot\boldsymbol{\rho}_{\mu}\right)\right.$$

$$\quad\left.-(m-g_{\sigma}\phi)\right]\Psi,$$ (2)

$$\mathcal{L}_{M} =\frac{1}{2}\left[\partial_{\mu}\sigma\partial^{\mu}\sigma-m_{\sigma}^{2}\sigma^{2}\right]$$

$$\quad-\frac{1}{4}\omega_{\mu\nu}\omega^{\mu\nu}+\frac{1}{2}m_{\omega}^{2}\omega_{\mu}\omega^{\mu}$$

$$\quad-\frac{1}{4}\boldsymbol{\varrho}_{\mu\nu}\cdot\boldsymbol{\varrho}^{\mu\nu}+\frac{1}{2}m_{\varrho}^{2}\boldsymbol{\varrho}_{\mu}\cdot\boldsymbol{\varrho}^{\mu},$$ (3)

$$\mathcal{L}_{NL} =-\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}$$

$$\quad+\Lambda_{\omega}g_{\varrho}^{2}\boldsymbol{\varrho}_{\mu}\cdot\boldsymbol{\varrho}^{\mu}g_{\omega}^{2}\omega_{\mu}\omega^{\mu},$$ (4)

where $\Psi$ represents the Dirac spinor nucleon doublet (proton and neutron) with a bare mass $m$, $g_{i}$ and $m_{i}$ are the coupling constants and the masses of the mesons $i=\sigma,\omega,\varrho$ and $\omega_{\mu\nu}=\partial_{\mu}\omega_{\nu}-\partial_{\nu}\omega_{\mu}$ and similar for $\boldsymbol{\varrho}_{\mu\nu}$. The equations of motion for $\sigma$, $\omega$ and $\varrho$ mesons are determined from the Euler-Lagrangian equations:

$$\sigma= \frac{g_{\sigma}}{m_{\sigma,\text{eff}}^{2}}\sum_{i}\rho_{i}^{S},$$ (5)

$$m_{\sigma,\text{eff}}^{2}=-m_{\sigma}^{2}-bmg_{\sigma}^{3}\sigma-cg_{\sigma}^{4}\sigma^{2},$$ (6)

$$\omega= \frac{g_{\omega}}{m_{\omega,\text{eff}}^{2}}\sum_{i}\rho_{i},$$ (7)

$$m_{\omega,\text{eff}}^{2}=m_{\omega}^{2}+\frac{\xi}{3!}g_{\omega}^{4}\omega^{2}+2\Lambda_{\omega}g_{\varrho}^{2}g_{\omega}^{2}\varrho^{2},$$ (8)

$$\varrho= \frac{g_{\varrho}}{m_{\varrho,\text{eff}}^{2}}\sum_{i}t_{3}\rho_{i},$$ (9)

$$m_{\varrho,\text{eff}}^{2}=m_{\varrho}^{2}+2\Lambda_{\omega}g_{\varrho}^{2}g_{\omega}^{2}\sigma^{2},$$ (10)

where $t_{3}=\pm 1/2$ is the isospin. The parameters $g_{\sigma}$, $g_{\omega}$, $\rho_{\varrho}$, $b$, $c$, $\xi$ and $\Lambda_{\omega}$ are sampled from Bayesian analysis.

To obtain the EOS in NSs conditions, we need to impose beta equilibrium and charge neutrality:

$$\mu_{p} =\mu_{n}-\mu_{e},$$ (11)

$$0 =\sum_{i=p,e,\mu}q_{i}\rho_{i},$$ (12)

where $q_{i}$ is the electric charge. However, this model should also satisfy the nuclear matter properties, constrained by Bayesian inference (see Sec. III). Symmetric nuclear matter (SNM) and pure neutron matter (PNM) equations are solved by imposing $\rho_{p}=\rho_{n}$ and $\rho_{p}=0$, respectively.

II.2 Quark phase

We used two different models to describe the quark phase: the Nambu-Jona-Lasinio (NJL) model [32, 33] and the Mean Field Theory of QCD (MFTQCD) [36]. The NJL model is a widely used model that includes all global QCD symmetries and reproduces chiral symmetry breaking in the vacuum. The second model is obtained by making approximations in the gluon field and exhibits behavior similar to that of the vectorial MIT bag model. For both models, we imposed chemical equilibrium and charge neutrality. A brief description of each model is provided below.

II.2.1 NJL

The Nambu-Jona-Lasinio (NJL) model [32, 33] is an effective model of point-like quark interactions. Despite the absence of gluons and a color confinement mechanism in this model, the NJL is well-suited for the description of large-density physics. This is due to its capacity to be designed to satisfy all the global symmetries of quantum chromodynamics (QCD) and to study manifestations of spontaneous chiral symmetry breaking [37]. In this work, the SU(3) NJL Lagrangian is given by

$$\mathcal{L} =\bar{\psi}(i\gamma_{\mu}\partial^{\mu}-m+\mu\gamma^{0})\psi$$

$$\quad+\frac{G}{2}\left[(\bar{\psi}\lambda_{a}\psi)^{2}+(\bar{\psi}i\gamma^{5}\lambda_{a}\psi)^{2}\right]$$

$$\quad+\kappa\left\{\det_{f}[\bar{\psi}(1+\gamma^{5})\psi]+\det_{f}[\bar{\psi}(1-\gamma^{5})\psi]\right\}$$

$$\quad+\mathcal{L}_{\text{int}},$$ (13)

where $m$ and $\mu$ are the quark current masses and chemical potential matrices, $\lambda_{a}$, with $a=1,2,\dots,8$, are the Gell-Mann matrices, and $\lambda_{0}$ is defined as $\lambda_{0}=\sqrt{2/3}\mathbb{1}$. The second term of Eq. 13 is the standard NJL term responsible for chiral symmetry breaking in the vacuum. The third one is implemented to explicitly break the $U(1)_{A}$ symmetry, as this is not a vacuum symmetry in QCD. The last term represents quark interaction terms added to better describe the physics of NSs. Here, we consider the following terms:

$$\mathcal{L}_{\text{int}} =-G_{\omega}\left[(\bar{\psi}\gamma^{\mu}\lambda_{0}\psi)^{2}+(\bar{\psi}\gamma^{\mu}\gamma^{5}\lambda_{0}\psi)^{2}\right]$$

$$\quad-G_{\rho}\sum_{a=1}^{8}\left[(\bar{\psi}\gamma^{\mu}\lambda_{a}\psi)^{2}+(\bar{\psi}\gamma^{\mu}\gamma^{5}\lambda_{a}\psi)^{2}\right]$$

$$\quad-G_{\omega\omega}\left[(\bar{\psi}\gamma^{\mu}\lambda_{0}\psi)^{2}+(\bar{\psi}\gamma^{\mu}\gamma^{5}\lambda_{0}\psi)^{2}\right]^{2}$$

$$\quad-G_{\sigma\omega}\sum_{a=0}^{8}\left[(\bar{\psi}\lambda_{a}\psi)^{2}+(\bar{\psi}i\gamma^{\mu}\gamma^{5}\lambda_{a}\psi)^{2}\right]$$

$$\quad\times\left[(\bar{\psi}\gamma^{\mu}\lambda_{0}\psi)^{2}+(\bar{\psi}\gamma^{\mu}\gamma^{5}\lambda_{0}\psi)^{2}\right]$$

$$\quad-G_{\rho\omega}\sum_{a=1}^{8}\left[(\bar{\psi}\gamma^{\mu}\lambda_{a}\psi)^{2}+(\bar{\psi}\gamma^{\mu}\gamma^{5}\lambda_{a}\psi)^{2}\right]$$

$$\quad\times\left[(\bar{\psi}\gamma^{\mu}\lambda_{0}\psi)^{2}+(\bar{\psi}\gamma^{\mu}\gamma^{5}\lambda_{0}\psi)^{2}\right].$$ (14)

To obtain the EOS, we apply the mean-field approximation. Effective mass and chemical potential are given by the gap equations

$$\tilde{m}_{i} =m_{i}+-2G\sigma_{i}-2\kappa\sigma_{j}\sigma_{k}$$

$$\quad+\frac{8}{3}G_{\sigma\omega}(\sigma_{i}^{2}+\sigma_{j}^{2}+\sigma_{k}^{2})\sigma_{i}$$ (15)

$$\tilde{\mu}_{i} =\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})$$

$$\quad-\frac{16}{9}G_{\omega\omega}(\rho_{i}+\rho_{j}+\rho_{k})^{3}$$

$$\quad-\frac{8}{3}G_{\sigma\omega}(\sigma_{i}^{2}+\sigma_{j}^{2}+\sigma_{k}^{2})(\rho_{i}+\rho_{j}+\rho_{k})$$

$$\quad-\frac{8}{9}G_{\rho\omega}(\rho_{i}+\rho_{j}+\rho_{k})$$

$$\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}),$$ (16)

with $i\neq j\neq k\in\{u,d,s\}$. For $T=0$, the grand canonical potential is given by

$$\Omega =\Omega_{0}+G\left(\sigma_{u}^{2}+\sigma_{d}^{2}+\sigma_{s}^{2}\right)+4\kappa\sigma_{u}\sigma_{d}\sigma_{s}$$

$$\quad-\frac{2}{3}G_{\omega}(\rho_{u}+\rho_{d}+\rho_{s})^{2}$$

$$\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})$$

$$\quad-\frac{4}{3}G_{\omega\omega}(\rho_{u}+\rho_{d}+\rho_{s})^{4}$$

$$\quad-4G_{\sigma\omega}\left(\sigma_{u}^{2}+\sigma_{d}^{2}+\sigma_{s}^{2}\right)\left(\rho_{u}+\rho_{d}+\rho_{s}\right)^{2}$$

$$\quad-\frac{8}{3}G_{\rho\omega}\left(\rho_{u}+\rho_{d}+\rho_{s}\right)^{2}$$

$$\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)$$

$$\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}},$$ (17)

where $\Omega_{0}$ is set to vanish the potential in the vacuum and $k_{F,f}=\sqrt{\tilde{\mu}_{f}^{2}+\tilde{m}_{f}^{2}}$ is the Fermi momentum. Here, we used the 3-momentum cutoff scheme ($\Lambda$). Imposing that the grand canonical potential must be stationary with respect to $\sigma_{i}$ and $\rho_{i}$ [39], i.e.,

$$\frac{\partial\Omega}{\partial\sigma_{i}}=\frac{\partial\Omega}{\partial\rho_{i}}=0,$$ (18)

we obtain

$$\sigma_{f} =-\frac{3}{\pi^{2}}\int_{k_{F_{f}}}^{\Lambda}dpp^{2}\frac{\tilde{m}_{f}}{E_{f}},$$ (19)

$$\rho_{f} =\frac{1}{\pi^{2}}k_{F_{f}}^{3},$$ (20)

at zero temperature.

The parameters $G$, $\kappa$, $m_{i}$ and $\Lambda$ are set to satisfy the mass and decay constant experimental data from $\pi^{\pm}$, $K^{\pm}$, $\eta$ and $\eta^{\prime}$ [49] (see Table 1). The coupling constants $G_{\omega}$, $G_{\rho}$, $G_{\omega\omega}$, $G_{\sigma\omega}$ and $G_{\rho\omega}$ are sampled from Bayesian analysis.

$\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 for NJL model, set to satisfy the experimental data from [49].

Furthermore, we set $P\to P+B$, where $B$ is a constant. This parameter exhibits behavior similar to the bag constant in the MIT bag model, strongly influencing the location of the phase transition point. This parameter is also sampled from Bayesian inference. The effects of the different multiquark interaction channels of Eq. (14) on the properties of hybrid stars, namely the interplay between the eight-quark vector interaction and the four-quark isovector-vector interaction, as well as higher-order repulsive interactions, have been studied in [50, 51, 52].

II.2.2 MFTQCD

In the mean field theory of quantum chromodynamics (MFTQCD), a decomposition of the gluon fields into low (soft) and high (hard) gluons is assumed in the QCD Lagrangian [36, 42, 43], i.e.,

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

where $G$ is the gluon field in momentum space and $A$ and $\alpha$ are the soft and hard gluon fields, respectively. Due to their small momenta, soft gluons are approximately constant and are replaced by their expected values in vacuum, given by [53, 54]

$$\left<A^{a\mu}A^{b\nu}\right> =-\frac{\delta^{ab}}{8}\frac{g^{\mu\nu}}{4}\mu_{0}^{2},$$ (22)

$$\left<A^{a\mu}A^{b\nu}A^{c\rho}A^{d\eta}\right> =\frac{\phi_{0}^{4}}{(32)(34)}\left[g_{\mu\nu}g^{\rho\eta}\delta^{ab}\delta^{cd}\right.$$

$$\quad\left.+g_{\mu}^{\rho}g_{\nu}^{\eta}\delta^{ac}\delta^{bd}+g_{\mu}^{\eta}g_{\nu}^{\rho}\delta^{ad}\delta^{bc}\right],$$ (23)

where $\phi_{0}$ and $\mu_{0}$ are energy scales to be determined. Assuming hard gluons have a large occupation number at all energy levels, they can be replaced by classical fields [55]

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

where $\alpha_{0}$ is a constant. The MFTQCD Lagrangian is obtained after a straightforward calculation

$$\mathcal{L}_{\text{MFTQCD}} =-B+\frac{m_{G}^{2}}{2}\alpha_{0}^{a}\alpha_{0}^{a}$$

$$\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},$$ (25)

where

$$m_{G}^{2} =\frac{9}{32}g^{2}\mu_{0}^{2},$$ (26)

$$B =\frac{9}{4(34)}g^{2}\phi_{0}^{4}=\left<\frac{1}{4}F^{a\mu\nu}F^{a}_{\mu\nu}\right>.$$ (27)

These are defined due to the fact that $m_{G}$ acts as a gluon mass in the MFTQCD Lagrangian and that $B$ exhibits behavior similar to the MIT bag constant. Using the energy-momentum tensor to calculate the equations of motion yields the EOS

$$P =\frac{27}{2}\xi^{2}\rho_{B}^{2}-B+P_{F},$$ (28)

$$\epsilon =\frac{27}{2}\xi^{2}\rho_{B}^{2}+B+\epsilon_{F},$$ (29)

where $P_{F}$ and $\epsilon_{F}$ are the pressure and energy density of a noninteracting Fermi gas of quarks and electrons, and $\xi\equiv g/m_{G}$. A more detailed deduction can be found in [36]. The final EOS exhibits behavior similar to that of a vectorial MIT bag model, in which the term $\xi$ acts as the vectorial term. In this work, the following values of masses were used: $m_{u}=5\text{ MeV}$, $m_{d}=7\text{ MeV}$ and $m_{s}=100\text{ MeV}$. The values of $\xi$ and $B$ are sampled by Bayesian inference.

III Bayesian Approach

Bayesian analysis samples the parameters of a model to satisfy the constraints which were imposed. This process is performed using Bayes’ theorem, given by

$$P(\text{model}|\text{data})=\frac{P(\text{data}|\text{model})P(\text{model})}{P(\text{data})},$$ (30)

where $P(\text{model}|\text{data})$ is the posterior distribution, $P(\text{data}|\text{model})$ the likelihood, $P(\text{model})$ the prior, and $P(\text{data})$ the evidence. In this work, we use the PyMultiNest [56, 57] sampler as part of the Bayesian Inference Library BILBY [58]. The sampler uses the nested sampling method, which starts with n random live points according to the prior distribution. In each iteration, the live point with the lowest likelihood is replaced by a new live point with a higher likelihood. This process moves the live points towards higher likelihood values. To apply this method, we must define the prior and likelihood distributions. The prior probability is the initial distribution of the parameters. Here, we use the uniform distribution defined in Table 2 (hybrid sets) and 3 (hadron set). We discuss five datasets: the hadronic set RMF, the hybrid sets NJL, NJL-GW, r-NJL and MFTQCD. The lowest and highest values of these distributions were chosen such that the posterior distributions did not show unjustified restrictions, except for the set r-NJL, for which the hadronic parameters were considered the same as the ones taken for the RMF set. This allows us to discuss, when building the hybrid EOS, the effect of forcing the hadronic prior space to coincide with the one considered for the hadronic EOS. Quark matter in sets NJL, NJL-GW, r-NJL is described by the NJL model and the quark model MFTQCD is used in the set MFTQCD.

To determine the likelihood - that is, the probability of obtaining the restrictions imposed for a particular parameter set - we use 10 different constraints. We divide them into 3 groups: experimental/observational data, guaranteeing a hybrid EOS, and corrections at M${}_{\text{max}}$.

Set NJL
NJL				RMF
Parameters	min	max		Parameters	min	max
$\xi_{\omega}$	0	0.5		$g_{\sigma}$	9	12
$\xi_{\rho}$	0	1		$g_{\omega}$	11	15
$\xi_{\omega\omega}$	0	30		$g_{\rho}$	9.546	15.000
$\xi_{\sigma\omega}$	0	8		BB	1.500	3.500
$\xi_{\rho\omega}$	0	50		CC	-4.627	-1.500
$B$ (MeV/fm3)	0	30		$\xi$	0	0.016
				$\Lambda_{\omega}$	0	0.103
Set MFTQCD
MFTQCD				RMF
Parameters	min	max		Parameters	min	max
$\xi_{Q}$ (MeV-1)	0	0.0018		$g_{\sigma}$	7	10
$B$ (MeV/fm3)	50	180		$g_{\omega}$	8	13
				$g_{\rho}$	8.000	15.000
				BB	1.000	9.000
				CC	-5.000	5.000
				$\xi$	0	0.040
				$\Lambda_{\omega}$	0	0.120
Set NJL-GW
NJL				RMF
Parameters	min	max		Parameters	min	max
$\xi_{\omega}$	0	0.5		$g_{\sigma}$	8	11
$\xi_{\rho}$	0	1		$g_{\omega}$	10	14
$\xi_{\omega\omega}$	0	30		$g_{\rho}$	9.546	15
$\xi_{\sigma\omega}$	0	8		BB	1.500	3.500
$\xi_{\rho\omega}$	0	50		CC	-4.627	-1.500
$B$ (MeV/fm3)	0	30		$\xi$	0	0.016
				$\Lambda_{\omega}$	0	0.103
Set r-NL
NJL				RMF
Parameters	min	max		Parameters	min	max
$\xi_{\omega}$	0	0.5		$g_{\sigma}$	8.010	9.691
$\xi_{\rho}$	0	1		$g_{\omega}$	9.084	12.167
$\xi_{\omega\omega}$	0	30		$g_{\rho}$	9.546	14.599
$\xi_{\sigma\omega}$	0	8		BB	2.205	6.903
$\xi_{\rho\omega}$	0	50		CC	-4.627	3.530
$B$ (MeV/fm3)	0	30		$\xi$	0	0.016
				$\Lambda_{\omega}$	0.036	0.103

Table 2: The lowest and highest values for the uniform distribution prior used for the data sets NJL, NJL-GW, MFTQCD and r-NJL. We defined $BB=b\times 10^{3}$ and $CC=c\times 10^{3}$.

Set RMF
Parameters	min	max
$g_{\sigma}$	6.5	13
$g_{\omega}$	6.5	15.5
$g_{\rho}$	6.5	16.5
BB	0.500	9.000
CC	-5.000	5.000
$\xi$	0	0.040
$\Lambda_{\omega}$	0	0.120

Table 3: The lowest and highest values for the uniform distribution prior used for the data set RMF. We defined $BB=b\times 10^{3}$ and $CC=c\times 10^{3}$.

1. Experimental/observational data

Nuclear matter properties ($\mathcal{L}_{\text{NMP}}$) constrain the EOS to satisfy the Nuclear matter properties (NMP) presented in Table 4 with the log-likelihood

$$\log(\mathcal{L}_{\text{NMP}})=-\frac{1}{2}\sum_{j}\left[\left(\frac{d_{j}-m_{j}(\theta)}{\sigma_{j}}\right)^{2}+\log(2\pi\sigma_{j}^{2})\right].$$ (31)

The first restriction in Table 4 concerns the saturation density ($\rho_{0}$). This is defined as the density of the symmetric nuclear matter at which the binding energy reaches its minimum, i.e.

$$\left.\frac{\partial(EA)}{\partial\rho}\right|_{\rho=\rho_{0}}=0,$$ (32)

where $EA=\epsilon/\rho-m_{n}$ with $m_{n}=939\text{ MeV}$. The second restriction comes from the known values of $EA(\rho=\rho_{0})$ and the incompressibility

$$K_{0}=9\rho_{0}^{2}\left.\frac{\partial^{2}(EA)}{\partial\rho^{2}}\right|_{\rho=\rho_{0}},$$ (33)

also for symmetric nuclear matter. The last one is applied to the symmetry energy at saturation,

$$J_{\text{sym,0}}=S(\rho_{0})=\frac{1}{2}\left.\frac{\partial^{2}(EA)}{\partial\delta^{2}}\right|_{\delta=0},$$ (34)

where $\delta=(\rho_{n}-\rho_{p})/\rho$ is the isospin asymmetry and the second derivative with respect to $\delta$ is taken at $\delta=0$, i.e. for symmetric matter .

It is also possible to calculate the skewness $Q_{0}$ ($n=3$) and kurtosis $Z_{0}$ ($n=4$) coefficients using

$$X_{0}^{(n)}=3^{n}\rho_{0}^{n}\left(\frac{\partial^{n}EA}{\partial\rho^{n}}\right)_{\delta=0},$$ (35)

and the symmetry energy slope $L_{\text{sym},0}$ ($n=1$), curvature $K_{\text{sym},0}$ ($n=2$), skewness $Q_{\text{sym},0}$ ($n=3$), and kurtosis $Z_{\text{sym},0}$ ($n=4$), given by

$$X_{\text{sym},0}^{(n)}=3^{n}\rho_{0}^{n}\left(\frac{\partial^{n}S(\rho)}{\partial\rho^{n}}\right)_{\rho_{0}}.$$ (36)

Pure Nuclear Matter ($\mathcal{L}_{\text{PNM}}$): constrains the EOS to satisfy the pure nuclear matter (PNM) energy per neutron from $\chi$EFT shown in Table 4 with the log-likelihood

$$\log(\mathcal{L}_{\text{PNM}})=\log\left[\prod_{j}\frac{1}{2\sigma_{j}}\frac{1}{\exp{\frac{|d_{j}-m_{j}(\theta)|-\sigma_{j}}{0.015}}+1}\right].$$ (37)

X-ray NICER Data ($\mathcal{L}_{\text{NICER}}$): The EOS must be able to describe observational data. Here we use the NICER data from J0030+0451 [59, 60] (with the ST+PDT hotspot model), J0740+6630 [61, 62] and J0437+4715 [63, 64]. The likelihood is given by

$$P(d_{X-ray}|EOS) =\int_{M_{\text{min}}}^{M_{\text{max}}}dmP(m|EOS)$$

$$\quad\times P(d_{X-ray}|m,R(m,EOS))$$

$$=\mathcal{L}_{\text{NICER}},$$ (38)

where

$$P(m|EOS)=\begin{cases}\frac{1}{M_{\text{max}}-M_{\text{min}}},\text{ if }M_{\text{min}}\leq m\leq M_{\text{max}},\\ 0,\text{otherwise},\end{cases}$$ (39)

with $M_{\text{min}}=1\text{M}_{\odot}$ and $M_{\text{max}}$ the maximum mass obtained from the EOS.

Gravitational Wave Data ($\mathcal{L}_{\text{GW}}$): We use the GW170817 data [65] from the LIGO-Virgo Collaboration. The likelihood is given by

$$\mathcal{L}_{\text{GW}}=\prod_{i}P(\Lambda_{1,i},\Lambda_{2,i},q_{i}|\mathcal{M}_{c},\boldsymbol{d}_{GW,i}(\boldsymbol{d}_{EM,i})),$$ (40)

where $\Lambda_{j,i}$ is the tidal deformability of the $j$ binary component, $q_{i}$ the mass ratio, $\mathcal{M}_{c}$ the chirp mass, and $\boldsymbol{d}_{GW,i}$ the observational data. The chirp mass is fixed at $\mathcal{M}_{c}=1.186\text{M}_{\odot}$.

NMP				PNM
Quantity	value/band	Ref.		Quantity	value/band	Refs.
$\frac{dEA}{d\rho}$	0			$EA_{PNM}(\rho=0.05)$	$6.8\pm 1.02$	[66]
$EA_{0}$	$-16\pm 0.02$	[67]		$EA_{PNM}(\rho=0.10)$	$10.5\pm 1.97$	[66]
$K_{0}$	$230\pm 30$	[68, 69]		$EA_{PNM}(\rho=0.15)$	$15.3\pm 3.44$	[66]
$J_{sym,0}$	$32.5\pm 1.8$	[70]

Table 4: Constraints for NMP and PNM.

2. Ensuring a hybrid EOS

Minimum distance between hadron and quark $P\times\mu_{B}$ ($\mathcal{L}_{\text{dist}}$): To generate a hybrid EOS through Maxwell construction, the quark and hadron equations should have an intersection point in the $P\times\mu_{B}$ plot. However, depending on the parameter values, there may be cases in which the $P_{Q}(\mu)$ and $P_{H}(\mu)$ curves do not intersect, and a phase transition is impossible. To ensure that we have this intersection point, we apply the following likelihood:

$$\log(\mathcal{L}_{\text{dist}})=-\frac{x^{2}}{2\delta^{2}},$$ (41)

where $x$ is the minimum distance between the $P_{Q}(\mu)$ and $P_{H}(\mu)$ curves and $\delta=0.01$. $\mathcal{L}_{\text{dist}}$ is a narrow Gaussian centered at $x=0$. The smaller the distance between the quark and hadron curves, the closer $\log(\mathcal{L}_{\text{dist}})$ is to zero.

Phase transition from hadron to quark ($\mathcal{L}_{\text{HtoQ}}$): In Maxwell’s construction, for each value of $\mu$, the pressure value will be the largest value between $P_{H}(\mu)$ and $P_{Q}(\mu)$. If $P_{H}(\mu)>P_{Q}(\mu)$ ($P_{Q}(\mu)>P_{H}(\mu)$), this indicates that the matter is in the hadron (quark) phase at $\mu$. A hybrid equation can thus be represented by the top plot of Fig. 1. Therefore, we will have hadron matter at $\mu<\mu_{\text{trans}}$, and quarks at $\mu>\mu_{\text{trans}}$. However, depending on the values of the parameters, we can have the case represented by the bottom plot of Fig. 1, resulting in a non-physical hybrid equation. To avoid this situation, we used the following likelihood:

$$\mathcal{L}_{\text{HtoQ}}=\frac{1}{1+\exp(ax+b)},$$ (42)

where $a=-6$, $b=1.5$, $x=P_{H}(\mu_{0})-P_{Q}(\mu_{0})$ with $\mu_{0}$ being a small value (in this case, $\mu_{0}=\mu_{Q}(\rho_{B}=0.235)$). This equation behaves as a smooth step function. Adjusting the $a$ and $b$ parameters, we can set the position of the step and its smoothness. Here, we choose $a$ and $b$ values to result in $\mathcal{L}_{\text{HtoQ}}\approx 1$ when $P_{H}(\mu_{0})>P_{Q}(\mu_{0})$ (hadron phase at $\mu_{0}$) and $\mathcal{L}_{\text{HtoQ}}\approx 0$ in the opposite case.

Figure 1: Example of a physical (top) and non-physical (bottom) phase transition. Cyan solid, black dashed and black dotted lines represent the hybrid, quark and hadron equations.

Range of phase transition ($\mathcal{L}_{\text{phT}}$): We constrained the density of the phase transition value using a super-gaussian centered at $\rho=0.275\,\text{fm}^{-3}$ with a standard deviation $\sigma=0.08$ and $p=5$. The likelihood is written as

$$\log(\mathcal{L}_{\text{phT}})=-\left[\frac{(\rho_{\text{trans}}-0.275)^{2}}{2(0.08)^{2}}\right]^{5}.$$ (43)

These chosen values imply that $\rho_{0}\lesssim\rho_{\text{trans}}\lesssim 0.40$.

3. Corrections at M${}_{\text{max}}$

Quarks inside M${}_{\text{max}}$ ($\mathcal{L}_{\text{Qmax}}$): In this study, we are interested in exploring cases in which quark matter can be found in stable NSs. To obtain these cases, the constraint $P_{\text{max}}>P_{\text{trans}}$ needs to be imposed. Fig. 2 shows the two possible hybrid mass-radius diagrams. The top plot shows the case in which we are interested, and the bottom plot shows the case we are avoiding. In the first case, NSs with $M>2\text{M}_{\odot}$ have quark matter in their core. In the second case, the maximum mass does not have enough pressure to deconfine matter. Hybrid stars are only possible in the unstable branch. To ensure that the EOS results in a maximum-mass NS with a significant amount of quark matter inside it, we used the smooth step function from Eq. (42) with $a=-0.2$, $b=20$, and $x=P_{\text{max}}-P_{\text{trans}}$. These values were chosen to position the step at approximately 100.

Figure 2: Example of a physical (top) and non-physical (bottom) phase transition. The red dot shows the maximum mass, and the green dot shows the phase transition.

Causality at M${}_{\text{max}}$ ($\mathcal{L}_{c_{s}^{2}\text{max}}$): To guarantee causality inside NSs, we impose that $c_{s}^{2}<1$ at the central density of $M_{\text{max}}$. To do this, we define the likelihood as a smooth step function, Eq. (42), with the values $a=100$, $b=-92$ and $x=c_{s}^{2}(\rho_{\text{max}})$.

Perturbative QCD ($\mathcal{L}_{\text{pQCD}}$): In [71], the authors developed a code [72] in which a Monte Carlo integration is performed for a given point of pressure, energy density, and baryonic density ($P,\epsilon,\rho_{B}$) to verify whether this point satisfies the pQCD constraint for an energy scale $X=[1/2,2]$. This code is based on the pQCD constraints discussed in [73]. The authors developed a method that constrains the $P\times\epsilon\times\rho_{B}$ space based solely on thermodynamic relations and ab initio calculations ($\chi$EFT [74] and pQCD [75]).

Full Likelihood

For the NJL and MFTQCD sets discussed in Sec. IV.1, and for the r-NJL set discussed in Sec. IV.3, the total likelihood can be written as

$$\log(\mathcal{L}) =\log(\mathcal{L}_{\text{NMP}})+\log(\mathcal{L}_{\text{PNM}})+\log(\mathcal{L}_{\text{NICER}})$$

$$\quad+\log(\mathcal{L}_{\text{dist}})+\log(\mathcal{L}_{\text{HtoQ}})+\log(\mathcal{L}_{\text{phT}})$$

$$\quad+\log(\mathcal{L}_{\text{Qmax}})+\log(\mathcal{L}_{c^{2}_{s}\text{max}})$$

$$\quad+\log[\mathcal{L}_{\text{pQCD}}(7\rho_{0})],$$ (44)

where we applied the pQCD constraint $\mathcal{L}_{\text{pQCD}}$ at $\rho_{B}=7\rho_{0}$. For the RMF set, also discussed in Sec. IV.1, the total likelihood is

$$\log(\mathcal{L}) =\log(\mathcal{L}_{\text{NMP}})+\log(\mathcal{L}_{\text{PNM}})+\log(\mathcal{L}_{\text{NICER}})$$

$$\quad+\log(\mathcal{L}_{c^{2}_{s}\text{max}})+\log[\mathcal{L}_{\text{pQCD}}(7\rho_{0})].$$ (45)

For the NJL-GW set discussed in Sec. IV.2, we applied the observation data GW170817 ($\mathcal{L}_{\text{GW}}$):

$$\log(\mathcal{L}) =\log(\mathcal{L}_{\text{NMP}})+\log(\mathcal{L}_{\text{PNM}})$$

$$\quad+\log(\mathcal{L}_{\text{NICER}})+\log(\mathcal{L}_{\text{GW}})$$

$$\quad+\log(\mathcal{L}_{\text{dist}})+\log(\mathcal{L}_{\text{HtoQ}})+\log(\mathcal{L}_{\text{phT}})$$

$$\quad+\log(\mathcal{L}_{\text{Qmax}})+\log[\mathcal{L}_{\text{pQCD}}(7\rho_{0})].$$ (46)

The value of $7\rho_{0}$ for the pQCD constraint was chosen due to the fact that the central densities at the maximum NS mass for the hybrid EOS sets can reach only $\sim 7\rho_{0}$ (see Table 5 and 6).

The evidence, which is only a normalization term, is estimated by the nested sampling method. We used 1,000 live points in all sets of equations and obtained 6514, 4879, 6037, 7521, and 5327 samples for the NJL, MFTQCD, RMF, NJL-GW and r-NJL sets, respectively.

IV Results

In this section, we present the EOSs, as well as other properties such as the mass-radius diagram, speed of sound, and trace anomaly, calculated using samples from Bayesian inferences. In total, we have obtained five EOS sets, four hybrid EOS sets, and one hadron EOS set. Three of the hybrid EOS sets use the NJL model (NJL all constraints except GW170817, NJL-GW also includes the GW170817 constraint discussed in Sec. IV.2, and r-NJL is determined with the same prior for the hadronic phase as the set RMF, not including GW170817), and the other uses the MFTQCD model for the quark phase. The same RMF model describes the hadron phase of all sets with, however, different priors (see Table 2). The hadronic priors have been chosen giving freedom to the Bayesian inference to search all the parameter space compatible with the constraints imposed through the likelihood, except for the r-NJL set. First, we compare three sets without the GW170817 constraint (NJL, MFTQCD and RMF). Next, we compare the two sets based on the NJL model (NJL and NJL-GW) to analyze the impact of the GW170817 constraint. In addition, we compare the NJL sets with two different hadron priors and without GW170817.

IV.1 Sets with different EOSs

Figure 3: Pressure versus energy density of the 90% of CI. Sets NJL (hybrid), MFTQCD (hybrid) and RMF (hadron) are represented in orange, cyan and hatched bands, respectively. Vertical lines indicate the 90% CI maximum for $\epsilon$. Band in gray represents the full (solid) and 90% of CI (dashed) of model-independent results from [76].

Sets NJL (hybrid), MFTQCD (hybrid), and RMF (hadron) are represented in Fig. 3 where the 90% credible intervals (CI) of pressure versus energy density for these three sets are plotted in orange, cyan, and hatched bands, respectively. When we compare the hybrid sets, we see that the MFTQCD set allows for a phase transition at low energy densities. This can also be seen in Table 5, which shows some numerical results, including the phase transition density. MFTQCD allows a phase transition at a density of $\sim 0.170\text{ fm}^{-3}$ (minimum 90% of CI), very close to the saturation density. In contrast, NJL allows a phase transition at a density of $\sim 0.304\text{ fm}^{-3}$ (minimum 90% of CI), approximately twice the saturation density. A similar result was obtained in [7], where the hadron phase was described by a fixed RMF equation, and Bayesian inference was applied only to the quark phase parameters. The three sets are compatible with the full model-independent results from [76], see the gray border in Fig. 3. However, the region defined in [76] was also constrained by the GW170817 detection, a condition that was not imposed on the NJL, MTFQCD, and RMF sets. In addition, the $\chi$EFT constraint considered in [76] is the NS matter pressure given in [74] which is different from the condition we apply in our analysis, the energy per neutron given in [66]. This explains why the NJL distribution may spread outside the range defined in [76]. Note that the RMF model follows the 90%CI band of [76], approximately. The MFTQCD model spans the 90%CI band and also covers a range below this band with the EOS already in the quark phase. The NJL model has a very hard hadronic EOS for densities after the $\chi$EFT band, spanning a region above the 90%CI band until the transition to the quark phase. The quark phase is compatible with the 90%CI band of [76]. Large quark cores are possible for particularly stiff hadronic EOSs.

Figure 4: Pressure (top) and speed of sound squared in units of $c^{2}$ (bottom) 90% CI distributions versus the baryonic density. Same color code from Fig. 3. Vertical bands represent the central density at the maximum NS mass of the 90% CI.

Fig. 4 shows the pressure and speed of sound squared versus baryonic density. The vertical bands represent the density at maximum NS mass. Interestingly, the speed of sound for each set is significantly different. For the NJL set, there are two bumps: the first is caused by the phase transition, and the second is caused by the appearance of the strange quark. After that, the speed of sound increases with density due to the term $\xi_{\omega\omega}$ [52, 7]. However, the EOSs are causal at the central density of maximum mass configurations (see the vertical band). The speed of sound squared takes values on the order of 0.5-0.6 for all models, at the center of maximum mass stars. The MFTQCD and RMF sets have similar speeds of sound at large densities, even though they have different components (quark and hadron respectively). They both increase until $\rho_{B}\approx 0.3\text{ fm}^{-3}$ (MFTQCD) and $\rho_{B}\approx 0.6\text{ fm}^{-3}$ (RMF), then stabilize at $c_{s}^{2}=0.5$. As shown by the vertical bands, the NJL set has smaller values for the central density ($\rho_{\text{max}}=0.947\text{ fm}^{-3}$ with maximum 90% CI), while the MFTQCD set can reach larger values for this quantity ($\rho_{\text{max}}=1.104\text{ fm}^{-3}$ with maximum 90% CI). See Table 5 for more numerical details.

Figure 5: Mass-radius diagram of the 90% of CI. Same color code from Fig. 3. Observational data are shown as PSR J0030+0451 (blue) [59, 60], PSR J0740+6620 (orange) [61, 62], PSR J0437+4715 (green) [63] and HESS J1731-347 (purple) [77] with $1\sigma$, $2\sigma$ and $3\sigma$ represented by solid, dashed and dotted lines, respectively.

The mass-radius diagram shown in Fig. 5 was obtained by solving the Tolman-Oppenheimer-Volkov (TOV) equations [78, 79]. The following observational data are represented in this figure with $1\sigma$ (solid), $2\sigma$ (dashed) and $3\sigma$ (dotted): PSR J0030+0451 (blue) [59, 60], PSR J0740+6620 (orange) [61, 62], and PSR J0437+4715 (green) [63, 64] by NICER, and HESS J1731-347 (purple) [77]. Interestingly, when we compare the hybrid sets with the hadron set, we see that the NJL set shifts the mass-radius diagram to the right, while the MFTQCD set shifts it to the left. Due to the small phase transition densities, the MFTQCD set can reach smaller radii than the other sets. It is the only set that is compatible with the HESS data (pink stain) at 68%. This result seems to indicate that the onset of quarks could occur at low densities. As a consequence, low mass NS could be hybrid stars. As we will discuss later, the lower the transition density to quark matter (and, therefore, the lower $M_{trans}$), the more compatible the MFTQCD curves are with the HESS data (see Fig. 18). As previously mentioned, the NJL model exhibits the opposite behavior, yielding larger radii. This occurs because the NJL EOS is not as soft as the MFTQCD EOS, and in order for the star to contain an appreciable amount of quark matter, the hadron phase should be described by a stiff EOS. To attain two solar mass stars, the term $\xi_{\omega\omega}$ plays an important role: at low densities its contribution is small, but its importance increases with density, making the quark EOS stiff enough. A stiff hadron EOS implies an increase in the radius of low mass stars. Note that the vector terms of the quark EOS allow for quite large radii at $\sim 2\text{M}_{\odot}$, compatible with the NICER results for J0740+6620. A more in-depth discussion can be found in the next section. The maximum masses of each set are $2.236\text{M}_{\odot}$, $2.315\text{M}_{\odot}$, and $2.185\text{M}_{\odot}$ at maximum 90% CI for the NJL, MFTQCD, and RMF sets, respectively. See Table 5 for more information. Notably, there is a increase in the maximum mass for the hybrid models. This is possible due to the vector terms present in the quark models. However, all sets can describe NICER data.

Figure 6: Corner plot of the hadron parameters for sets NJL, MFTQCD and RMF. Same color code from Fig. 3.

Fig. 6 shows the corner plot of the hadron parameters. It can be seen that the MFTQCD model does not significantly change the values of the hadron parameters compared to the RMF set, except for the $\xi$ parameter, to which the model is not sensitive. This is because this term plays a role at high densities when the phase transition to quark matter has already occurred. Within the NJL model, a different result was obtained: the values of $g_{\sigma}$, $g_{\omega}$ are much larger, and there is no superposition with the values obtained for RMF and MFTQCD; $g_{\rho}$ shows a wider distribution and may take larger values than those of RMF and MFTQCD; BB, CC and $\Lambda$ peak at smaller values. $\xi$ is the only parameter for which its posterior remains approximately the same. The differences in the parameter distributions explain why the NJL set has larger radii. It is a consequence of imposing that two solar mass stars are described and that the hybrid star has a non-negligible quark core. To attain these conditions, the hadron EOS must be quite stiff to allow an early transition to quark matter. However, one consequence is that the NJL set does not satisfy the GW170817 constraints, as can be seen in Fig. 7. This is not the case with the MFTQCD and RMF sets, which satisfies the GW170817 constraints, although they have not been included in the Bayesian inference. Fig. 7 shows the relation between the binary mass ratio $q=M_{2}/M_{1}<1$ and the effective tidal deformability, given by

$$\tilde{\Lambda}=\frac{16}{13}\frac{(12q+1)\Lambda_{1}+(12+q)q^{4}\Lambda_{2}}{(1+q)^{5})},$$ (47)

where $\Lambda_{i}$ and $M_{i}$ are the tidal deformability and mass of the i-th NS in the binary system, respectively. For the different sets, $\tilde{\Lambda}$ values were calculated for $0.73<q<1$ and a fixed value of $M_{\text{chirp}}=(M_{1}M_{2})^{3/5}/(M_{1}+M_{2})^{1/5}=1.186\text{M}_{\odot}$, in accordance with event GW170817 [65]. The observational data from the LIGO/Virgo collaboration for event GW170817 is shown in the figure by the solid green contours (50% CI), dashed contours (90% CI) and dotted contours (99% CI) [65]. The NJL set does not satisfy the GW170817 constraint due to the large radii it predicts. In the next section, we will apply the GW170817 constraint and analyze the differences caused by this restriction. One direct consequence of the new constraint is to restrict the maximum radius in this set.

Figure 7: 90% of CI of tidal deformability. Same color code from Fig. 3. The GW170817 observational data from the LIGO/Virgo collaboration is shown by solid (50% CI), dashed (90% CI) and dotted (99% CI) green contours [65].

Figure 8: 90% of CI of measure of conformability ($d_{c}$), polytropic index ($\gamma$) and trace anomaly ($\Delta$) versus baryonic density. Same color code from Fig. 3.

Fig. 8 shows the polytropic index $\gamma=d\ln P/d\ln\epsilon$, the trace anomaly $\Delta=1/3-P/\epsilon$, and the measure of conformability $d_{c}=\sqrt{\Delta^{2}+\Delta^{\prime 2}}$, with $\Delta^{\prime}=d\Delta/d\ln\epsilon$, as functions of the baryonic density. These quantities were used in model-independent EOSs to identify quark matter inside NSs [80], [81], [82]. In [81], it was proposed that $\gamma<1.75$ could indicate deconfined phase behavior. In [80], $d_{c}<0.2$ was used instead. However, none of our sets follow these trends. In $d_{c}$ plots, we can identify the phase transition by the bump at $\rho\approx 0.2-0.4\text{ fm}^{-3}$ for both hybrid sets. After the phase transition, the value of $d_{c}$ decreases and only reaches $0.2$ at densities of approximately $\rho\approx 0.5\text{ fm}^{-3}$ and $\rho\approx 0.7\text{ fm}^{-3}$ for the NJL and MFTQCD sets, respectively. Similar behavior is observed in the $\gamma$ plots. Furthermore, the hadron set reaches the $d_{c}$ and $\gamma$ limits at $\rho\approx 0.8\text{ fm}^{-3}$, even though it contains no quarks. This may suggest that $d_{c}$ and $\gamma$ are not suitable quantities for indicating the presence of quark matter by themselves. All sets reach negative values of $\Delta$ at large densities.

	Set NJL			Set NJL-GW
quant	median	min	max	median	min	max
$\rho_{\text{trans}}$	0.353	0.304	0.388	0.362	0.314	0.391
$P_{\text{trans}}$	57.466	35.143	78.623	55.989	33.814	77.024
$\epsilon_{\text{H,trans}}$	350.549	295.526	391.007	358.957	304.433	393.36
$\epsilon_{\text{Q,trans}}$	395.697	331.009	468.817	395.183	334.473	458.06
$\Delta\epsilon_{\text{trans}}$	44.971	22.492	91.833	36.207	19.691	76.780
$M_{\text{max}}$	2.130	2.018	2.236	2.108	1.993	2.212
$R_{\text{max}}$	12.466	11.716	13.122	12.147	11.526	12.827
$\rho_{\text{max}}$	0.829	0.733	0.947	0.871	0.767	0.968
$\epsilon_{\text{max}}$	995.402	857.809	1178.842	1053.422	898.290	1211.631
$c_{s,\text{max}}^{2}$	0.529	0.329	0.776	0.564	0.340	0.795
$M_{\text{Q,max}}$	1.103	0.757	1.511	1.176	0.777	1.555
$R_{\text{Q,max}}$	7.955	6.956	8.933	8.069	7.026	8.952
$R_{1.4\text{M}_{\odot}}$	13.695	13.390	14.051	13.448	13.193	13.698
$M_{\text{trans}}$	1.728	1.341	2.008	1.634	1.253	1.940
	Set MFTQCD			Set r-NJL
quant	median	min	max	median	min	max
$\rho_{\text{trans}}$	0.222	0.170	0.308	0.369	0.333	0.395
$P_{\text{trans}}$	5.752	2.451	15.851	54.239	36.589	70.712
$\epsilon_{\text{H,trans}}$	208.795	157.329	295.494	364.841	324.201	395.105
$\epsilon_{\text{Q,trans}}$	289.862	254.155	337.472	399.421	349.543	446.832
$\Delta\epsilon_{\text{trans}}$	77.373	24.691	122.923	33.315	20.009	60.217
$M_{\text{max}}$	2.133	1.970	2.315	1.996	1.863	2.154
$R_{\text{max}}$	11.273	10.685	11.954	12.199	11.547	12.847
$\rho_{\text{max}}$	0.976	0.852	1.104	0.861	0.757	0.962
$\epsilon_{\text{max}}$	1205.841	1052.149	1361.615	1013.698	863.204	1191.993
$c_{s,\text{max}}^{2}$	0.487	0.458	0.515	0.388	0.200	0.744
$M_{\text{Q,max}}$	2.018	1.781	2.224	1.010	0.572	1.432
$R_{\text{Q,max}}$	10.266	9.405	11.101	7.818	6.471	8.633
$R_{1.4\text{M}_{\odot}}$	12.111	11.571	12.649	13.265	13.050	13.436
$M_{\text{trans}}$	0.327	0.174	0.630	1.547	1.248	1.771

Table 5: The 90% CI of the following quantities: the phase transition density ($\rho_{\text{trans}}$, in fm^-3), pressure ($P_{\text{trans}}$, in MeV/fm³), hadronic energy density ($\epsilon_{\text{H,trans}}$, in MeV/fm³) and quark energy ensity ($\epsilon_{\text{Q,trans}}$, in MeV/fm³); measure of the strength of the phase transition ($\Delta\epsilon_{\text{trans}}=\epsilon_{\text{Q,trans}}-\epsilon_{\text{H,trans}}$, in MeV/fm³); maximum mass ($M_{\text{max}}$, in M_⊙) and radius ($R_{\text{max}}$, in km), density ($\rho_{\text{max}}$, in fm^-3), energy density ($\epsilon_{\text{max}}$, in MeV/fm³), speed of sound ($c_{s,\text{max}}^{2}$) quark core mass ($M_{\text{Q,max}}$, in M_⊙) and quark core radius ($R_{\text{Q,max}}$, in km) at the maximum NS mass; radius of the $1.4\text{M}_{\odot}$ ($R_{1.4\text{M}_{\odot}}$, in km); mass of stars with central pressure equal to $P_{\text{trans}}$ ($M_{\text{trans}}$, in M_⊙). Results of the hybrid sets.

Set RMF
quant	median	min	max
$M_{\text{max}}$	2.039	1.905	2.185
$R_{\text{max}}$	10.761	10.339	11.195
$\rho_{\text{max}}$	1.103	1.003	1.212
$\epsilon_{\text{max}}$	1354.565	1191.639	1533.405
$c_{s,\text{max}}^{2}$	0.582	0.478	0.692
$R_{1.4\text{M}_{\odot}}$	12.297	11.923	12.714

Table 6: The 90% CI of the following quantities: maximum mass ($M_{\text{max}}$, in M_⊙) and radius ($R_{\text{max}}$, in km), density ($\rho_{\text{max}}$, in fm^-3) and energy density ($\epsilon_{\text{max}}$, in MeV/fm³), speed of sound ($c_{s,\text{max}}^{2}$) at the maximum NS mass; radius of the $1.4\text{M}_{\odot}$ ($R_{1.4\text{M}_{\odot}}$, in km). Results of the hadronic set.

IV.2 Set with GW170817 constraint

In the previous section, we have identified discrepancies in the hadron phase parameter values and the radius on the mass-radius diagram between the NJL set and the other sets. In this section, we compare the NJL set with the NJL-GW set, which is built from the same hadron and quark models (NJL model for the quark phase and RMF model for the hadron phase) and constraints as the NJL set, but with the addition of the GW170817 constraint (see Sec. III).

Figure 9: Pressure versus energy density of the 90% of CI. NJL, NL-GW, r-NJL and RMF sets are represented in orange, dotted purple, blue and hatched bands, respectively. Vertical lines indicate the 90% CI maximum for $\epsilon$. Solid and dashed gray contours represent the full and 90% Cl result from [76].

Fig. 9 shows the pressure as a function of energy density. The orange band represents the 90% CI of the NJL set. The purple band with a dotted pattern represents the 90% CI of the NJL-GW set. This figure also shows the results for the RMF set (hatched pattern) and the r-NJL set (blue), which will be discussed in Sec. IV.3. In the hadron phase region ($\epsilon\lesssim 400\text{ MeV/fm}^{3}$), the NJL-GW set exhibits lower pressure than NJL and falls inside the [76] envelope. However, in the quark phase region ($\epsilon\gtrsim 400\text{ MeV/fm}^{3}$), the two sets have almost completely overlapping 90% CIs. Fig. 10 shows the corner plot for the quark phase parameters. All of the quark parameters have approximately the same posterior distribution. $\xi_{\omega\omega}$, an important parameter for obtaining $M_{\text{max}}>2\text{M}_{\odot}$ [7, 51, 52], has slightly smaller values for the NJL-GW.

Figure 10: Corner plot of the quark phase parameters for the NJL (orange) and the NJL-GW (purple) sets.

Fig. 11 shows the speed of sound in terms of baryonic density. The orange, purple dotted, hatched and blue bands represent the results for the NJL, NJL-GW, r-NJL and RMF sets, respectively. The r-NJL will be discussed in Sec. IV.3. Vertical bands show the central density of maximum-mass NSs. Notably, NJL, NJL-GW and r-NJL sets exhibit a double bump in the speed of sound, as discussed before: the first representing the quark-hadron phase transition and the second representing the appearance of strange quarks. The phase transition bump occurs at slightly lower densities for the NJL set than for the NJL-GW set because the NJL hadronc phase is stiffer and therefore the phase transition density to quark matter is lower for NJL than for NJL-GW (see Table 5). However, this difference is less noticeable for the second bump, and results essentially due to the small differences of the couplings associated to the flavor dependent terms. The NJL set reaches higher speeds of sound due to its larger values of $\xi_{\omega\omega}$.

Figure 11: Speed of sound versus density of the 90% of CI. Vertical bands represent the central density at the maximum NS mass. Same color code from Fig. 9.

Fig. 12 shows the 90% CI of the mass-radius diagram. Interestingly, even though these hybrid sets use the same models for the hadron and quark phases, the mass-radius relations differ significantly. The NJL-GW set has a notably smaller radius than the NJL set. The most significant difference is the maximum 90% CI of $R_{1.4\text{M}_{\odot}}$. The difference between the maximums is 0.353 km, while the difference between the minimums is $\sim$0.197 km (see Table 5). This implies that the constraint imposed by GW170817 restricts the radius from increasing too much. This difference is explained by the correlation between $g_{\omega}$ and $R_{1.4\text{M}_{\odot}}$, as seen in Fig. 13. Notably, $R_{1.4\text{M}_{\odot}}$ increases with $g_{\omega}$. This occurs because a larger $g_{\omega}$ gives rise to a stiffer hadron EOS. For the NJL-GW set, $g_{\omega}$ cannot reach as large values as for the NJL set. However, the minimum value of this coupling remains almost the same. The coupling $g_{\sigma}$ also changes because it is linearly correlated with $g_{\omega}$. This correlation is imposed by the binding energy at saturation.

Figure 12: Mass-radius diagram of the 90% of CI. Same color codes from Fig. 9 and same observational data from Fig. 5.

Figure 13: Corner plot of $g_{\sigma}$, $g_{\omega}$ and $R_{1.4}$ for the NJL sets with different priors.

Even when applying the GW170817 constraint, it is notable that NJL-GW’s radius is still larger than those of the MFTQCD and RMF sets. At this point, we might ask, “Why do Bayesian analyzes prefer EOS with larger radii for the NJL sets?” Larger values of $g_{\omega}$, which increase the radius, decrease the value of $\xi_{\omega}$ (see Fig. 10) and, therefore, that of the transition density $\rho_{\text{trans}}$. This allows for larger values of the $\xi_{\omega\omega}$ coupling constant, which increases the maximum mass. In other words, increasing $g_{\omega}$ offsets the increase in $\xi_{\omega}$, resulting in values of $\rho_{\text{trans}}$ within the Bayesian constraint. Additionally, current observational data set stronger constraints on the mass than on the radius, implying that Bayesian inference will prioritize larger $M_{\text{max}}$ than $R_{1.4}\sim 13$ km.

Figure 14: Corner plot of $\xi_{\omega\omega}$, $g_{\sigma}$, $g_{\omega}$, $\rho_{\text{trans}}$ and $M_{\text{max}}$ for the NJL sets with different priors.

Fig. 15 shows the relation between the binary mass ratio ($q$) and the effective tidal deformability ($\tilde{\Lambda}$) for the NJL, NJL-GW, r-NJL and RMF sets. Even though the GW170817 constraint was included in the Bayesian analysis, the NJL-GW set still struggles to describe the observational data. Nevertheless, the results for NJL-GW are within the 99% CI of the data and are more favorable than those of the NJL set. The NJL equation has difficulty reconciling the tidal deformability data, which requires a softer equation, with the two solar mass data, which requires a stiffer equation.

Figure 15: 90% of CI probability distributions of the NS of the tidal deformability. Same color codes from Fig. 9. The GW170817 observational data from the LIGO/Virgo collaboration is shown by solid (50% CI), dashed (90% CI) and dotted (99% CI) green contours [65].

IV.3 Hadron prior

In this subsection, we briefly discuss the choice of the hadron couplings prior and compare NJL, NJL-GW and r-NJL sets. NJL and NJL-GW sets have been introduced and discussed. The r-NJL set was obtained by considering the same constraints as those for the NJL set, but with a more restricted prior for the hadron coupling. For this set, we used a uniform distribution with a minimum and maximum of the 90% CI from Set 0 of [31]. In that paper, Bayesian inference was applied to the same RMF model used in this work, which constrained the NMP and PNM. In this subsection, we compare a prior that does not restrict parameter values, as in the NJL and NJL-GW sets, with a prior obtained from the previous Bayesian analysis done in [31] that did not consider quarks.

Fig. 9 shows the 90% CI of pressure times the energy density for the NJL (orange), NJL-GW (purple dotted), r-NJL (blue) and RMF (hatched) sets. When comparing the hybrid sets, the r-NJL has the lowest pressure in the hadron phase ($\epsilon\lesssim 400\text{MeV/fm}^{3}$). However, it is notable that this set is stiffer than the RMF set. For large energy densities, the r-NJL set is the softest — it has the lowest pressure compared to the other three sets, including the RMF set.

The speed of sound as a function of density is shown in Fig. 11. The r-NJL phase transition bump occurs at larger densities than the NJL and NJL-GW sets due to its larger value of $\rho_{\text{trans}}$ (see Table 5). At higher densities, it is interesting to note that the r-NJL has a broader range and can reach lower values than the RMF model. However, the r-NJL’s 90% maximum is smaller than those of the NJL and NJL-GW.

Fig. 12 shows the mass-radius diagram. The r-NJL set describes the smallest radii compared to the other NJL hybrid sets. However, it also has the smallest mass of all sets, with $M_{\text{max}}^{\text{med}}=1.996\text{M}_{\odot}$. Although the r-NJL set is mostly stiffer than the RMF set (see Fig. 9), the fact that the RMF set is stiffer at high densities is sufficient for it to reach higher mass values.

In Fig. 15, already discussed above, we also include the result for r-NJL. It is interesting to note that the restricted prior of the hadron parameters (r-NJL), has a greater effect on the $\tilde{\Lambda}$ values than by applying the GW170817 restriction (NJL-GW). Although the results of the r-NJL set only describe the GW170817 event at 99% CI, it is the set that uses the NJL equation that best satisfies this observational data. However, it is also the set with the lowest maximum mass values, as we can see from Tables 5 and 6.

V Implications

In this section, we compare some of the NS and NMP properties predicted by the different models. We first discuss the NS properties.

Distributions of the $\rho_{\text{trans}}$ (top plot) and $\Delta\epsilon_{\text{trans}}$ (lower plot) are shown in Fig. 16. Although the constraint on the transition density $\rho_{\text{trans}}$ allows transitions between $0.16\lesssim\rho_{\text{trans}}\lesssim 0.40\text{ fm}^{-3}$, we notice that sets using the NJL model allow phase transitions only above $\sim 0.25\text{ fm}^{-3}$ (all NJL sets). However, the MFTQCD set can describe a phase transition within the entire allowed range. Interestingly, the MFTQCD set prefers a lower phase transition with a median of $\rho_{\text{trans}}=0.222\text{ fm}^{-3}$, while the NJL sets prefer a higher phase transition density, with a median of $\rho_{\text{trans}}=0.353\text{ fm}^{-3}$, $0.362\text{ fm}^{-3}$ and $0.369\text{ fm}^{-3}$, respectively, for the NJL, the NJL-GW and the r-NJL sets. The NJL-GW distribution is only slightly shifted to larger densities when compared with NJL, r-NJL showing a larger shift. The behavior of the NJL models with respect to that of the MFTQCD indicates that the NJL EOSs are stiffer, favoring a late phase transition.

The $\Delta\epsilon_{\text{trans}}$ is defined as $\Delta\epsilon_{\text{trans}}=\epsilon_{\text{Q, trans}}-\epsilon_{\text{H,trans}}$, where $\Delta\epsilon_{\text{trans}}$ is a measure of the strength of the phase transition. The $\Delta\epsilon_{\text{trans}}$ distribution has the widest distribution for the MFTQCD set and peaks at larger $\Delta\epsilon_{\text{trans}}$, indicating that the strongest phase transitions occur with the MFTQCD model for quark matter. NJL, NJL-GW and r-NJL sets have peaks for quite smaller $\Delta\epsilon_{\text{trans}}$, $\sim 30$ MeV/fm³, about three times smaller than the corresponding $\Delta\epsilon_{\text{trans}}$ for MFTQCD.

Figure 16: Histogram of $\rho_{\text{trans}}$ and $\Delta\epsilon_{\text{trans}}$ for the NJL, MFTQCD and NJL-GW sets.

Figure 17: Histogram of $\rho_{\text{trans}}$ and $M_{\text{max}}$ for the NJL, MFTQCD, RMF, NJL-GW and r-NJL sets. The maximum M${}_{\text{max}}$ reached by each set is $2.384M_{\odot}$, $2.444M_{\odot}$, $2.302M_{\odot}$, $2.327M_{\odot}$, $2.291M_{\odot}$, respectively.

Fig. 17 shows the histogram of the maximum NS mass (top), the mass of the star for which quarks start to nucleate, i.e. with $P_{c}=P_{\text{trans}}$ (middle) and the quark core mass of the maximum mass NS (lower). Analyzing the maximum mass values, the NJL, MFTQCD, RMF reach similar results, with a median of $M_{\text{max}}=2.130,2.133,2.039\text{M}_{\odot}$, respectively. MFTQCD is reaching the largest masses, 2.444 $\text{M}_{\odot}$, while NJL-GW, r-NJL and RMF only reach $\sim 2.384,\,2.327$ and 2.291 $\text{M}_{\odot}$. It is interesting to see that the hybrid stars reach the largest masses. This is due to the fact that while the properties of the hadronic EOS are strongly constrained by the NMP, the quark EOS has more freedom, being constrained by causality and pQCD constraints. Although the NJL and MFTQCD sets have almost the same median, the MFTQCD set has a slightly wider distribution. The smaller maximum mass attained by NJL-GW with a median of $M_{\text{max}}=2.108\text{M}_{\odot}$, is due to the smaller values of $\xi_{\omega\omega}$ of this model, as discussed in section IV.2. See Table 5 for more numerical results.

The middle plot of Fig. 17 shows the mass of stars with a central pressure equal to $P_{\text{trans}}$ ($M_{\text{trans}}$), i.e., the lowest mass at which quark matter can be found. Within the MFTQCD set, quark matter can exist in NSs with very low masses, below 1 $\text{M}_{\odot}$, the median value being equal to 0.327 $\text{M}_{\odot}$. These are essentially quark stars with a hadronic crust. NJL models predict the presence of quarks inside heavier stars with a mass above 1 $\text{M}_{\odot}$. For NJL, NJL-GW and r-NJL sets, the median mass is 1.728 $\text{M}_{\odot}$, 1.634 $\text{M}_{\odot}$ and 1.547 $\text{M}_{\odot}$, respectively, all above 1.5 $\text{M}_{\odot}$. Identifying such different compositions, one expects that other NS properties, such as those obtained from cooling, NS modes or binary neutron star mergers, will distinguish the two scenarios. Notice that in Fig. 3, the pressure range covered by the MFTQCD model below the 90% CI band corresponds to these low mass stars with a quark core. See Table 5 for more information.

Figure 18: Speed of sound as a function of density (left), mass-radius diagram (middle) and the binary mass ratio versus the effective tidal deformability (right) for the NJL-GW (top), NJL (middle) and MFTQCD (bottom) sets. 90% CI of equations with $M_{\text{trans}}<0.5\text{M}_{\odot}$, $0.5\text{M}_{\odot}<M_{\text{trans}}<1.4\text{M}_{\odot}$ and $M_{\text{trans}}>1.4\text{M}_{\odot}$ are represented by the blue, dotted and red bands, respectively.

Therefore, the hybrid sets suggest that PSR J0740+6620 could be a hybrid star and that PSR J0030+0451 is inconclusive. Interestingly, the NJL-GW set has larger values for $\rho_{\text{trans}}$ but smaller values for $M_{\text{trans}}$ compared to the NJL set. This is because the mass of a NS is obtained from the TOV equations, which describe the variation of pressure (rather than density) in terms of radius. Comparing the $P_{\text{trans}}$ values in Table 5 reveals that the NJL set has larger $P_{\text{trans}}$ values than the NJL-GW set.

The quark core mass of $M_{\text{max}}$ is shown in the lower plot of Fig. 17. Due to the small values of $\rho_{\text{trans}}$, MFTQCD set has the heaviest quark core, allowing NSs with a quark core mass larger than $2\text{M}_{\odot}$. These NSs are made mostly of quark matter. The NJL, NJL-GW and r-NJL have similar distributions, with a median of $M_{\text{Q,max}}^{\text{med}}=1.103$, $1.176\text{M}_{\odot}$ and $1.010\text{M}_{\odot}$, respectively. This would be about half of the maximum NS mass.

	NJL-GW	NJL	MFTQCD
$M_{\text{trans}}<0.5\;\text{M}_{\odot}$	0	0	4113
$0.5<M_{\text{trans}}<1.4\;\text{M}_{\odot}$	1112	536	766
$M_{\text{trans}}>1.4\;\text{M}_{\odot}$	6409	5978	0
total	7521	6514	4879

Table 7: Number of equations for each band in Fig. 18.

Figure 19: Histogram of nuclear matter properties for the NJL, MFTQCD, RMF, NJL-GW and r-NJL sets.

Fig. 18 shows the influence of $M_{\text{trans}}$ on the speed of sound (left), the mass-radius diagram (middle), and the effective tidal deformability (right). The three hybrid sets - NJL-GW (top), NJL (middle), and MFTQCD (bottom) - were separated in the cases of $M_{\text{trans}}<0.5\text{M}_{\odot}$, $0.5\text{M}_{\odot}<M_{\text{trans}}<1.4\text{M}_{\odot}$ and $M_{\text{trans}}>1.4\text{M}_{\odot}$. The 90% CI of these cases are represented by the blue, dotted, and red bands, respectively. Due to the lack of equations with $M_{\text{trans}}<0.5\text{M}_{\odot}$ for the NJL-GW and NJL, and $M_{\text{trans}}>1.4\text{M}_{\odot}$ for the MFTQCD sets, these cases are omitted. Table 7 shows the number of equations for each band: most of the MFTQCD (NJL and NJL-GW) stars have $M_{\text{trans}}<0.5\,\text{M}_{\odot}$ ($M_{\text{trans}}>1.4\,\text{M}_{\odot}$).

For the NJL-GW and NJL sets, the speed of sound for $0.5\text{M}_{\odot}<M_{\text{trans}}<1.4\text{M}_{\odot}$ has a smaller first bump than that of the $M_{\text{trans}}>1.4\text{M}_{\odot}$ case due to a phase transition at a lower density. The second bump caused by the onset of strangeness is not affected by the value of $M_{\text{trans}}$. In the mass-radius diagram, smaller $M_{\text{trans}}$ values correspond to smaller radii, and therefore, for the MFTQCD model the $M_{\text{trans}}<0.5\text{M}_{\odot}$ set can describe the HESS data. The same behavior is seen with tidal deformability, implying that for the models discussed, lower phase transitions are more compatible with GW170817 data.

Set RMF
quant	median	min	max
$\rho_{0}$	0.159	0.158	0.161
EA	-16.000	-16.032	-15.968
$K_{0}$	250.940	228.158	283.916
$Q_{0}$	-445.834	-510.475	-346.521
$J_{\text{sym,0}}$	32.138	29.678	34.717
$L_{\text{sym,0}}$	46.898	34.300	66.925
$K_{\text{sym,0}}$	-140.665	-184.537	-69.975
$Q_{\text{sym,0}}$	1136.710	407.590	1563.186

Table 8: 90% CI of NMP of the hadronic RMF set.

	Set NJL			Set NJL-GW
quant	median	min	max	median	min	max
$\rho_{0}$	0.159	0.158	0.161	0.159	0.158	0.161
EA	-16.000	-16.033	-15.967	-16.000	-16.032	-15.968
$K_{0}$	229.623	188.253	269.753	228.829	185.733	271.410
$Q_{0}$	417.230	-85.649	1663.222	71.412	-273.811	573.546
$J_{\text{sym,0}}$	31.829	29.363	34.470	31.149	28.706	33.708
$L_{\text{sym,0}}$	62.443	52.668	74.218	54.256	48.117	62.705
$K_{\text{sym,0}}$	110.178	17.644	206.472	91.128	6.448	169.655
$Q_{\text{sym,0}}$	902.673	-426.280	1476.560	1002.040	202.071	1506.580
	Set MFTQCD			Set r-NJL
quant	median	min	max	median	min	max
$\rho_{0}$	0.159	0.158	0.161	0.159	0.158	0.160
EA	-16.000	-16.034	-15.967	-16.000	-16.032	-15.968
$K_{0}$	240.731	218.267	266.119	248.612	213.777	278.984
$Q_{0}$	-557.802	-619.236	-462.206	-97.085	-307.695	73.622
$J_{\text{sym,0}}$	32.563	30.206	35.048	31.473	28.863	33.655
$L_{\text{sym,0}}$	49.115	36.964	69.491	51.564	44.450	59.747
$K_{\text{sym,0}}$	-158.530	-191.779	-106.057	2.801	-49.118	78.403
$Q_{\text{sym,0}}$	976.226	337.799	1475.125	1276.510	906.791	1557.227

Table 9: 90% CI of NMP of the hybrid sets: NJL, NJL-GW, MFTQCD and r-NJL.

$M$	$1.2\text{M}_{\odot}$		$1.4\text{M}_{\odot}$		$1.6\text{M}_{\odot}$		$1.8\text{M}_{\odot}$
$dM/dR$	+	-	+	-	+	-	+	-
NJL	6306	208	6008	506	4896	1618	2509	4005
NJL-GW	7268	253	6656	865	4699	2822	1653	5868
r-NJL	5180	147	4562	765	2358	2969	206	5121
MFTQCD	4392	487	4326	553	3459	1420	1119	3760
RMF	822	5215	175	5862	70	5967	38	5999

Table 10: Number of EOSs with positive/negative slope in the mass-radius diagram at $M=1.2,1.4,1.6,1.8\text{M}_{\odot}$.

Fig. 19 shows the NMP distributions for NJL, NJL-GW, r-NJL, MFTQCD and RMF sets. Tables 9 and 8 show the 90% CI of NMP of the hybrid and hadronic sets, respectively. In the different Bayesian analyzes considered, the values of $\rho_{0}$, EA, $K_{0}$ and $J_{\text{sym},0}$ are constrained, implying that their distributions are similar. While the other quantities shown in Fig. 19 are not constrained, it is notable that MFTQCD and RMF generally have similar distributions, with a tendency of MFTQCD to have slightly smaller properties. This might be expected because the hadron phase properties of MFTQCD are defined for densities of about 0.15-0.25 fm^-3, which is close to the density that defines NMP. Additionally, the stiffness required to describe a two solar mass system is defined by the quark phase.

The NJL-GW and r-NJL sets show NMP distributions closer to those of the MFTQCD and RMF sets than the NJL set. It is interesting to note that the MFTQCD and RMF sets prefer negative values for $K_{\text{sym},0}$, which is in agreement with [83]. The r-NJL takes values between RMF and NJL, NJL-GW and $K_{\text{sym},0}$ is essentially centered at zero. The NJL and NJL-GW sets take positive values of $K_{\text{sym},0}$ in agreement with predictions of [84], where a RMF model that also includes the $\delta$-meson has been fitted to the CREX and PREX data. In this paper, positive values of $K_{\text{sym},0}$ were proposed to reconcile the PREX and CREX results.

The symmetry energy slope $L_{\text{sym},0}$ controls both the neutron skin thickness and the radius of low-mass nuclei. In [85], a value of $L_{\text{sym},0}=58.7\pm 28.1\text{MeV}$ was obtained at 1$\sigma$ from both experimental and observational data. All sets are in agreement with this result.

We next analyze the slope of the mass-radius curve at some given NS masses. In [86], it has been shown that the behavior of this quantity has a direct influence on the maximum mass attained and on the NS radius. In addition, it has been discussed in [87], where the slopes of the mass-radius curves from nucleonic and hyperonic EOSs were analyzed, that the composition could affect the sign of the slope. Table 10 shows the number of EOSs with positive and negative slope, i.e. $dM/dR$, at $M=1.2,1.4,1.6,1.8\text{M}_{\odot}$. The RMF set is the set that presents the largest number of negative slope values for all masses: 86%, 97%, 99% and 99.5% for stars with $M=1.2,1.4,1.6$ and 1.8 $\text{M}_{\odot}$. The hybrid sets mostly have positive values for $M=1.2,1.4,1.6$ $\text{M}_{\odot}$. In particular, for $M=1.2\,\text{M}_{\odot}$, about 97% of the EOSs in all NJL sets have a positive slope, and 90% in the MFTQCD EOSs. These values decrease to $\sim 90\%$ for the NJL models and $M=1.4\,\text{M}_{\odot}$, while it remains almost unchanged for MFTQCD. For $M=1.8\,\text{M}_{\odot}$ there are still about 20%-38% EOSs with a positive slope. Only r-NJL shows a smaller value, $\sim 3\%$ in line with the RMF set. In [87] similar conclusions have been drawn for the nucleonic EOS: essentially, a negative slope for all masses was obtained. For the hyperonic EOSs, positive values were mostly obtained for small masses, similarly to the results obtained for the hybrid sets in this work. Within both studies, a positive value of $dM/dR$ at small masses could indicate the presence of exotic degrees of freedom, such as quarks or hyperons. A significant difference in our present study is that at $1.8\text{M}_{\odot}$ there is still a reasonable fraction of curves with a positive slope. This reflects the two model description of the hybrid stars, with the quark EOSs being less constrained. The vector terms are responsible for the stiff EOSs that justify the positive slope for this large mass. If in the future it is confirmed that two solar mass stars like J0740–6620 (with $R=12.48^{+1.28}_{-0.88}$ at 68% CI [61]) have a larger radius than 1.4 solar mass stars as J0614–3329 (with $R=10.29^{+1.01}_{-0.86}$ at 68% CI [88]), a quark core of strongly interacting matter could explain this radius difference.

Figure 20: Compactness $C=M_{\text{max}}/R_{\text{max}}$ distribution for all sets: NJL-GW (purple), NJL (orange), MFTQCD (cyan), r-NJL (blue) and RMF (black dashed). Upper limit $C=1/3$ [89] is shown as the vertical red line.

In [89], the authors discuss the maximum compactness of NSs. Considering an agnostic description of the EOS and imposing both observational and theoretical constraints, they conclude that if the maximum compactness $\mathcal{C}_{max}$ is attained by the maximum mass configuration, then $\mathcal{C}_{\text{max}}=1/3$. In our study, we consider EOS motivated by microscopic models and impose similar constraints . It is natural that the extremes are not attained because the models are not general enough, however, it is still interesting to identify the maximum values obtained for the compactness, defined as $C=M_{\text{max}}/R_{\text{max}}$, within our five datasets and analyze the properties of the corresponding stars. Fig. 20 shows its distribution for all five generated sets. The r-NJL model has the lowest compactness values. This is due to its small maximum mass combined with a large radius. The NJL and NJL-GW models have the second and third smallest values of compactness, respectively. They are greater than the r-NJL because these sets have greater maximum masses. However, they are lower than the MFTQCD and RMF due to their large radii. The MFTQCD and RMF sets have the largest compactness values, reaching 0.290 and 0.299, respectively (see Table 11). All sets satisfy the $C\leq 1/3$ constraint, represented by the vertical red line, obtained in [89].

Set	$C_{\text{min}}$ (M/R)	$C_{\text{max}}$ (M/R)
NJL-GW	0.219 (1.93/13.00)	0.291 (2.29/11.46)
NJL	0.219 (1.96/13.17)	0.288 (2.21/11.31)
MFTQCD	0.263 (1.86/10.41)	0.290 (2.44/12.39)
r-NJL	0.209 (1.87/13.18)	0.286 (2.27/11.72)
RMF	0.246 (1.87/11.17)	0.299 (2.26/11.17)

Table 11: Minimum and maximum values of compactness for each set. In parentheses, we have the respective values of mass and radius in M_⊙ and km.

A corner plot of $C$, $M_{\text{max}}$, and $R_{\text{max}}$ is given in Fig. 21. The $C\times M_{\text{max}}$ panel shows that there is a correlation between these two quantities: larger values of $M_{\text{max}}$ result in larger values of $C$. This behavior is expected, given the definition of $C=M_{\text{max}}/R_{\text{max}}$. The opposite behavior is expected for the correlation between $C$ and $R_{\text{max}}$, as seen with the NJL, NJL-GW and r-NJL sets. However, RMF does not have this correlation. For the MFTQCD set, $C$ increases with $R_{\text{max}}$. This is because these two sets have a correlation between $M_{\text{max}}$ and $R_{\text{max}}$, differing from the NJL, NJL-GW and r-NJL sets.

Figure 21: Corner plot of the compactness ($C$), maximum mass (${M_{\text{max}}}$) and radius at the maximum mass ($R_{\text{max}}$).

Fig. 22 shows the mass-radius diagram of the maximum (solid) and minimum (dashed) of $C=M_{\text{max}}/R_{\text{max}}$ for each set. The mass-radius curves of the maximum (minimum) $C$ values have large (small) $M_{\text{max}}$, respectively. This is a behavior also found in [89]. Two different trends are obtained for the configurations with the minimum compactness. NJL models show the lowest values mainly because these models predict large radii. In our approach the maximum NS mass is defined by the pulsar J0740-6620 and we get values below 2 $M_{\odot}$. Taking as minimum $M_{TOV}=2.2M_{\odot}$, as the authors in [89], we would have obtained $\sim 0.24$. For the other two data sets the minimum compactness is attained for stars with a small radius and the value of the compactness does not drop below 2.46 for RMF and 2.63 for MFTQCD. These results are also compatible with the discussion in [89] where the minimum compactness was obtained for large radii. In Fig. 23, the plot shows the distribution of the maximum compactness as a function of the maximum mass. All points fall above the lower limit obtained in [89]

Figure 22: Mass-radius diagram of the maximum (solid lines) and minimum (dashed lines) of $C=M_{\text{max}}/R_{\text{max}}$.

Figure 23: Compactness as function of $M_{\text{max}}$. The red dashed line represents the lower limit of $C$ obtained in [89].

VI Conclusions

This work explores the possibility of quark matter existing inside NSs. We used microscopic models to derive the EOSs. The relativistic mean field (RMF) model was used to describe the hadron phase of matter. For the quark phase, we employed the Nambu-Jona-Lasinio (NJL) model and the mean field theory of quantum chromodynamics (MFTQCD). We applied a Maxwell construction to build the quark-hadron phase transition. We obtained large sets of EOSs using Bayesian inference with theoretical constraints from nuclear matter properties and pQCD calculations, as well as observational constraints from NICER. In order to obtain hybrid stars with a large quark core, a minimal pressure difference between the pressure at the phase transition and that at the maximum NS mass configuration is imposed. Five sets of EOSs were obtained in total: four with hybrid EOSs and one with only hadron EOSs (the RMF set). Three of the hybrid sets used the NJL model for the quark phase: two of the sets, NJL and NJL-GW, differed because the second one was also constrained by the GW170817 observation. In addition, we have built a third NJL set to discuss the effect of the hadron phase prior. The other set, the MFTQCD set, used the MFTQCD model to describe the quark core.

We have considered microscopic models to describe hybrid stars; as a result, the composition of these objects is known, and it is possible to discuss the effect of the presence of a quark core. The hadron phase is described by a framework that has been quite successful in describing the properties of finite nuclei and nuclear matter. A larger uncertainty concerns the description of quark matter; therefore, two different approaches were considered, which resulted in different compositions. In addition, no exotic quark phases, such as color-superconducting phases, have been considered. The treatment of the quark phase and the phase transition requires further discussion. It is, however, interesting to identify the regions in the pressure–energy density plane spanned by the different sets and compared to the expectations from agnostic approaches, as discussed in [91, 76]. While the RMF EOS distribution follows the 90% CI predicted in [76], the hadron phase of the NJL sets lies above the 90% CI although the quark phase is consistent with the 90% CI of the agnostic study. Concerning the MFTQCD, the EOS probabilities of the low mass hybrid stars predicted by this model fall partially in a region below the 90% CI of the agnostic study.

When comparing the different hybrid models, we conclude that the MFTQCD set allows the phase transition to occur at low densities with $\rho_{\text{trans}}^{\text{min}}=0.170$ fm^-3 (90% CI), essentially giving rise to quark stars with a crust, whereas the NJL, NJL-GW and r-NJL sets only allow the phase transition to occur above $2\rho_{0}$ with $\rho_{\text{trans}}^{\text{min}}=0.304{\text{ fm}^{-3}}$, $\rho_{\text{trans}}^{\text{min}}=0.314{\text{ fm}^{-3}}$ and $\rho_{\text{trans}}^{\text{min}}=0.333{\text{ fm}^{-3}}$, respectively (90% CI). Note that similar onset densities were obtained in [51, 52], where the DDME2 was considered to describe the nucleonic EOS and the NJL, with the terms we are considering for the quark phase. We have imposed a constraint on the phase transition density value, allowing it to occur in $0.16\lesssim\rho_{\text{trans}}\lesssim 0.40\text{ fm}^{-3}$. The low $\rho_{\text{trans}}$ values imply a large quark core mass. MFTQCD set can describe NSs with a quark core mass of $M_{\text{Q,max}}^{\text{med}}=2.018\text{M}_{\odot}$, while the NJL, NJL-GW and r-NJL, all with phase transitions above $2\rho_{0}$, predict a quark core mass with $M_{\text{Q,max}}^{\text{med}}=1.103\text{M}_{\odot}$, $M_{\text{Q,max}}^{\text{med}}=1.176\text{M}_{\odot}$ and $M_{\text{Q,max}}^{\text{med}}=1.010\text{M}_{\odot}$, respectively.

Concerning the maximum mass, the RMF set has the smallest $M_{\text{max}}$, reaching $({M_{\text{max}}})^{\text{max}}=2.185\text{M}_{\odot}$ (90% CI). The presence of quarks results in slightly larger maximum masses. The MFTQCD, NJL and NJL-GW sets can reach $({M_{\text{max}}})^{\text{max}}=2.315\text{M}_{\odot}$, $({M_{\text{max}}})^{\text{max}}=2.236\text{M}_{\odot}$ and $({M_{\text{max}}})^{\text{max}}=2.212\text{M}_{\odot}$ (90% CI), respectively. This occurs because there is some freedom in building the quark phase, which is only constrained by causality and the pQCD constraints, together with a maximum mass of at least two solar masses. Both quark models considered contain vector contributions that can stiffen the EOS inside NSs while still satisfying the pQCD constraints at very large densities. This stiffening could explain the larger radius predicted for the two solar mass pulsar J0740–6620 compared with the radius of the 1.4 solar mass pulsar J0614–3329.

We also analyzed the possibility of quark matter in a $1.4\text{M}_{\odot}$ NS. Due to the small phase transition densities, the MFTQCD set indicates that NS with $1.4\text{M}_{\odot}$ have quark matter in their inner core. This set results in $({M_{\text{trans}}})^{\text{max}}=0.630\text{M}_{\odot}$ (90% CI). However, the NJL and NJL-GW sets favor EOSs without quark matter inside $1.4\text{M}_{\odot}$ NS. While these sets do not exclude this possibility, as they obtain $({M_{\text{trans}}})^{\text{min}}=1.341\text{M}_{\odot}$ and $({M_{\text{trans}}})^{\text{min}}=1.253\text{M}_{\odot}$ (90% CI), the median values indicate that quark matter is present only for $M>1.728\text{M}_{\odot}$ and $M>1.634\text{M}_{\odot}$, respectively. Although the presence of quark matter in $1.4\text{M}_{\odot}$ NSs is inconclusive, all hybrid sets agree that quark matter is present in NSs with masses greater than 2 solar masses.

The slope of the mass-radius curves may also carry some information about the possible existence of exotic matter inside NSs. It was shown that if the slope of the mass-curve is negative for all masses, there is a large probability that the NS does not contain exotic matter, while the opposite conclusion is drawn with a positive slope at 1.2 or 1.4 $\text{M}_{\odot}$.

Within our sets, we could also conclude that a polytropic index below 1.75 or the trace anomaly related quantity $d_{c}<0.2$ does not necessarily indicate the presence of quark matter, as suggested in [81, 80]: some of our hybrid EOS predict $d_{c}>0.2$ in the quark phase, and the hadronic EOS may have $d_{c}<0.2$ at large densities. We have also analyzed the maximum NS compactness determined from our datasets. Values below 0.3 were obtained, consistent with the findings of [89].

All sets can describe both theoretical constraints and observational data. However, NS data carry large uncertainties, especially regarding the radius values. This makes it difficult to make strong statements about the matter phase of their inner cores. Third-generation telescopes are expected to measure radii within $<100$ meters [92], which would impose strong constraints on the EOS. This could provide more information about the composition of matter and its properties under the extreme conditions that are uniquely reproduced by NSs.