PNJL model at zero temperature: three-flavor case

We propose a three-flavor version of the Polyakov-Nambu-Jona-Lasino (PNJL) model at zero temperature regime, by implementing a traced Polyakov loop ($\Phi$) dependence in the scalar, vector and 't Hooft channel strengths. We study the thermodynamics of this model, named as PNJL0, with special attention for the first order confinement/deconfinement phase transition for which $\Phi$ is the order parameter. For the symmetric quark matter case, an interesting feature observed is a strong reduction of the constituent strange quark mass ($M_s$) at the chemical potential related to point where deconfinement takes place. The emergence of $\Phi$ favors the restoration of chiral symmetry even for the strange quark. We also investigate the charge neutral system of quarks and leptons in weak equilibrium. As an application, we construct a hadron-quark phase transition with a density dependent hadronic model coupled to the SU(3) PNJL0 model. In this case, the quark side is composed by deconfined particles. This approach is used to determine mass-radius profiles compatible with recent data from the Neutron Star Interior Composition Explorer (NICER) mission.

We propose a three-flavor version of the Polyakov-Nambu-Jona-Lasino (PNJL) model at zero temperature regime, by implementing a traced Polyakov loop (Φ) dependence in the scalar, vector and 't Hooft channel strengths. We study the thermodynamics of this model, named as PNJL0, with special attention for the first order confinement/deconfinement phase transition for which Φ is the order parameter. For the symmetric quark matter case, an interesting feature observed is a strong reduction of the constituent strange quark mass (Ms) at the chemical potential related to point where deconfinement takes place. The emergence of Φ favors the restoration of chiral symmetry even for the strange quark. We also investigate the charge neutral system of quarks and leptons in weak equilibrium. As an application, we construct a hadron-quark phase transition with a density dependent hadronic model coupled to the SU(3) PNJL0 model. In this case, the quark side is composed by deconfined particles. This approach is used to determine mass-radius profiles compatible with recent data from the Neutron Star Interior Composition Explorer (NICER) mission.

I. INTRODUCTION
The understanding of the physics related to particles interacting strongly is a fundamental issue studied for both, theoretical and experimental researchers. Many aspects of the strongly interacting matter can be directly or indirectly accessed from experiments, as for example, the charge-parity violation (CPV) observed in heavy mesons decays. The reader is addressed to Ref. [6] and references therein for an overview including the role of the strong force in forming the observed CPV pattern.
From the theoretical point of view, on the other hand, the physics of quarks and gluons is described by Quantum Chromodynamics (QCD) [7][8][9]. Its nonperturbative regime, i.e., QCD at large distances, or equivalently, low energies, can be investigated through the lattice discretization of QCD in the Euclidean space [10][11][12] where the numerical calculations are based on Monte Carlo simulations [13]. However, the so called fermion sign problem [14] restricts such technique when quarks at finite chemical potential are included in the system. A second approach applied to study the infrared region of QCD is the use of continuous methods relying on Dyson-Schwinger equations (DSE) in Euclidean space [15], used to determine the equations of motion of the n-point functions. Another practical method of treating quarks and gluons, in systems in which such particles are the funda-mental degrees of freedom, is the use of phenomenological models. Probably, the simplest one used for this aim is the famous MIT bag model [16]. In its original version, it establishes that quarks are inside a bag that effectively corresponds to a confining potential. The constant B is added to the free quark pressure, and all thermodynamics is based on this simple assumption.
Other widely used effective quark model is the Nambu-Jona-Lasinio (NJL) one [17][18][19][20][21]. It was proposed by Y. Nambu and G. Jona-Lasinio before QCD and originally described nucleons and pions. In the context of the quark interaction, it presents some similarities with QCD such as the dynamical breaking of chiral symmetry, and the associated quark mass generation due to the quark condensates. Similar nonperturbative phenomena also occur in the Bardeen-Cooper-Schrieffer (BCS) model of superconductivity [22] (in which the NJL model was inspired), where interactions give rise to energy gaps. In the NJL model, the "gap" in the quark mass is produced by the point-like interactions between Dirac particles.
Although capable of describing some important features of QCD, the NJL model fails to reproduce a particular phenomenon, namely, the infrared dynamics of confinement as well as the deconfinement of quarks, related to the so called "asymptotic freedom", discovered by Gross, Wilczek and Politzer [23]. In the NJL model there is no information regarding the interaction among the gluons. Such an issue was corrected by K. Fukushima [24], who proposed the inclusion of these mediators through a background field obtained from the trace of a matrix in color space associated to the gluon gauge field A 4 . The so-called traced Polyakov loop, Φ (and its conjugate Φ * ), mimics the gluon dynamics in an effective way. Physically, Φ can be seen as the order parameter of the confinement/deconfinement phase transition, since it is related to e −F/T , with F being the free energy of an external static quark. In the confined phase, the free energy of a single quark is infinite and, consequently, Φ = 0. In the deconfined phase, on the other hand, F is finite and Φ = 0. However, we remind that Φ is not enough to have a true confinement in the PNJL model at finite temperature since it only presents a statistical confinement while asymptotic quark states still exists. The aforementioned analysis originally does not apply at zero temperature in the PNJL model and in most of its parametrizations, since in this case the Polyakov potential vanishes and the modified Fermi-Dirac distributions reduces to the step functions. Because of that, confinement effects are lost and it is not possible to investigate high density deconfined quark systems at T = 0. In order to circumvent this issue, it was proposed in Refs. [25,26] a modification of the strengths of the couplings in the NJL model by making them dependent on Φ. As a direct consequence, the Polyakov loop potential became nonvanishing even at T = 0. In this way, the new model constructed for flavor SU(2) quarks, named as PNJL0 model, was able to exhibit the confinement/deconfinement phenomenology at T = 0. In this work, we extend these previous studies by generalizing the PNJL0 model to include flavor SU(3) quarks. We investigate its thermodynamics at T = 0 and also perform an application in the context of the hadron-quark phase transition.
The paper is organized as follows. In Sec. II we review the main equations of the traditional PNJL at finite temperature. In Sec. III, we generalize the PNJL0 model at T = 0 by taking into account the dynamics of the strange s quark. The thermodynamics of the SU(3) model is studied for both cases, namely, symmetric quark matter, and charge neutral system of quarks and leptons in weak equilibrium. In the former, we verify a strong reduction of the constituent strange quark mass (M s ) at the chemical potential where the finite value of Φ emerges (deconfined phase). An application to hybrid stars is also performed. The results indicate that it is possible to generate hybrid stars compatible with the recent data from the Neutron Star Interior Composition Explorer (NICER) mission. Finally, summary and concluding remarks are presented in Sec. IV.

II. THREE-FLAVOR POLYAKOV-NAMBU-JONA-LASINIO MODEL
We start by presenting the Lagrangian density of the three-flavor version of the PNJL model, that includes a scalar and a vector 4-fermion interaction. It reads (qγ µ λ a q) 2 + (qγ µ γ 5 λ a q) 2 where q is a vector of three spinors q f for f = u, d, s, m = diag(m u , m d , m s ) is a matrix of current quark masses in flavor space, and λ a are the SU(3) Gell-Mann matrices. The last term, a determinant in flavor space, corresponds to the so called 't Hooft interaction [27] and is responsible for the large value of the η ′ mass due to the anomaly U A (1) [28]. Such a kind of interactions were firstly considered by Kobayashi and Maskawa [29]. The differences between this model and the NJL one [17][18][19][20][21][30][31][32] are the replacement of ∂ µ by D µ ≡ ∂ µ + iA µ , where A µ = δ µ 0 A 0 and A 0 = gA 0 a λ a /2 (g is the gauge coupling), and the inclusion of the Polyakov potential U(Φ, Φ * , T ). It depends on the traced Polyakov loop and its conjugate, Φ and Φ * , with Φ is defined in terms of and DS10 [37,38]. The traced Polyakov loop is found from ∂Ω PNJL /∂Φ = 0, and the current quark masses M u,d,s are determined through solutions of Eqs. (7), in which condensates ρ s u,d,s given in Eqs. (5) are also used. Finally, from Eq. (4) all the other thermodynamical quantities can be determined, namely, pressure, entropy and energy densities. The expressions are obtained from P PNJL = −Ω PNJL , S PNJL = −∂Ω PNJL /∂T , and E PNJL = T S PNJL − P PNJL + f µ f ρ f , respectively.

A. Formulation
At this point, we remind the reader that the the PNJL model equations, previously presented, are reduced to those related to the NJL one at zero temperature regime. The origin of this result is twofold. First, one has that the generalized Fermi-Dirac distributions, given in Eqs. (8) and (9), become step functions θ(k F f − k) at T = 0 (k F f is the Fermi momentum of the quark f ). Therefore, Φ disappears in all momentum integrals. The second reason is due to the gluonic contribution enclosed by the Polyakov potential. Notice that U(Φ, 0) vanishes for the most known versions. In this case, Eq. (4) reads In order to circumvent this problem, we implemented, in Ref. [26], the introduction of Φ in the SU(2) NJL model at T = 0 by requiring that the strengths of the scalar and vector channels are vanishing in the deconfined phase, i.e., at Φ = 1. In that approach we incorporate expected effects from QCD that at low densities, and corresponding large inter-particle distances, quarks should interact strongly while at short distance the interaction should be weakened. The former is associated to the nonperturbative infrared physics from QCD that enhances the interaction between the effective degrees of freedom, related to the quarks in the model. The latter should represent the ultraviolet physics of perturbative QCD where the quarks interact weakly due to the asymptotic freedom phenomenon. Therefore, this physics is incorporated in a dynamical way such that at low densities quarks interact strongly and at large densities weakly towards to the deconfinement phase transition. For this purpose, we use the traced Polyakov loop as an effective scalar background field at the zero temperature regime. Although this quantity is not strictly defined at T = 0, we use it to bring the complexity of QCD dynamics by incorporating, in a effective way, the transition from strong (infrared) to weak (ultraviolet) regimes of the quark-quark interaction. Here, we extended this phenomenology also to the three-flavor version of the model. The modifications in the couplings for this case are and Such changes lead to with By choosing to incorporate the terms G s Φ 2 , G V Φ 2 and KΦ 2 of the grand canonical potential into U defined in Eq. (22), it is possible to identify Eq. (21) as the zero temperature version of Eq. (4). As in Ref. [26], we refer to this construction as the PNJL0 model, now written in its SU(3) version. Quark masses and chemical potentials are given, respectively, by and where and Notice that in the three-flavor version of the PNJL0 model, the modification of the couplings, Eqs. (18)- (20), also induces the definition of a new Polyakov potential, Eq. (22), in which the backreaction of quarks in the gluonic sector happens, as well as the inverse one, namely, gluons affecting quarks. This last effect is intrinsic in the original PNJL model but the former is absent. In the SU(3) PNJL0 model, the backreaction is complete, i.e., each sector interacts each other, exactly as in its twoflavor version [26]. We remark that the modifications pointed out in Eqs. (18)- (20) lead to the Polyakov po- However, this potential is not able to generate nonvanishing solutions for Φ when the condition ∂Ω PNJL0 /∂Φ = 0 is applied. Therefore, we add the term U 0 given in Eq. (16) in the definition of U and verify that it is responsible to ensure Φ = 0. Moreover, it is also important in order to restrict the traced Polyakov loop in the range of 0 Φ 1 [37,38]. Finally, the pressure and energy density of the SU(3) PNJL0 model, obtained from Eq. (21), are written as and respectively. The constant Ω vac = −P vac is added to the pressure and energy density in order to ensure that in vacuum (ρ u = ρ d = ρ s = 0) one has Ω PNJL0 (ρ f = 0) = −P PNJL0 (ρ f = 0) = 0 and E PNJL0 (ρ f = 0) = 0.

B. Thermodynamics of the model (symmetric quark matter)
We remark that from now on, all results are obtained at zero temperature regime.
Since the main equations of the model are defined, we are able to explore some thermodynamical features of the SU(3) PNJL0 model. Firstly, we need to define values for its free parameters, namely, G s , G V , K, m u , m d , m s and Λ. Here we use the RKH parametrization [19,47] in which G s = 3.67/Λ 2 , K = −12.36/Λ 5 , m u = m d = 5.5 MeV, m s = 140.7 MeV, and Λ = 602.3 MeV. The last quantity is the cutoff parameter. As in the NJL model, the PNJL0 one is nonrenormalizable. Therefore, it is needed to adopt a certain regularization scheme in order to avoid divergent contributions in the momentum integrals. Here we use a three-momentum cutoff for this purpose [19,48,49]. This parameters set produces the values of f π = 92.4 MeV, m π = 135 MeV, m K = 497.7 MeV, m η = 514.8 MeV, and m ′ η = 957.8 MeV, respectively, for the pion decay constant and mass, and for the masses of pseudoscalar mesons K, η, and η ′ [19,47]. For the vacuum, it also leads to M (vac) With regard to the strength of the vector channel, we remind that if a Fierz transformation of the color-current interaction is performed, the value of G V = 0.5G s is found [19,20]. In principle, this "canonical" value should be used in all calculations involving NJL/PNJL models with the vector channel included. However, G V is often taken as a free parameter. Some authors use this freedom to fix the vector meson masses, such as the ρ meson [50]. Other ones also point out to changes in the magnitude of G V due to medium effects. For instance, instanton-anti-instanton pairing effects at high temperature can change its magnitude [51]. Furthermore, it is also known that G V plays a significant hole in the phase diagram of the strongly interacting matter since it shifts the location of the critical end point (CEP) [35,52]. The increasing of G V induces the disappearance of the CEP and, consequently, the vanishing of the first order phase transition line. Therefore, as in Ref. [26], we adopt here the approach of setting G V as a free parameter in order to verify its influence in our calculations.
Finally, the gluonic sector of the model present two more constants, namely, T 0 and a 3 , appearing in the last term of the Polyakov potential, Eq. (22). The former is fixed to 190 MeV [26,53]. Concerning the constant a 3 , we remind that it appears in the quantity denoted by U 0 , Eq. (16), firstly used in Ref. [37]. In those works, the authors remark that the free parameters presented in the Polyakov Potential U, that contains U 0 , are determined in order to reproduce lattice data, as well as known information about the phase diagram. In this case, the value obtained by them is a 3 = −0.4. However, it is worth to notice that their model is completely different from the PNJL0 model considered here. If we take a 3 = −0.4, no solutions of Φ = 0 are found, as we will show in the first figure later on. Therefore, one has G V and a 3 as free parameters of the SU(3) PNJL0 model, as well as they are in its two-flavor version [26].
At zero temperature, one can also use the thermodynamical relation given by in order to write the quark chemical potentials as i.e., in terms of the baryonic, isospin, and strangeness chemical potentials, namely, µ B , µ I , and µ S , respec-tively. Firstly, we work at the so called symmetric quark matter, the one in which all quark potentials are the same, namely, situation found by taking µ I = µ S = 0 in Eq. (30). At this point, we remark to the reader that the concept of symmetric matter in this context is different from the one used in nuclear matter described by hadronic models [54], in which nucleons are in equal number and present identical chemical potentials, since they have the same rest masses and densities. In the case of the three-flavor quark system studied here, m u = m d but m s is quite different. Therefore, densities of the light quarks are equal each other but different from the strange quark one, despite the validity of Eq. (31).
As a first investigation of the symmetric quark matter described by the SU(3) PNJL0 model, we show in Fig. 1 Ω PNJL0 as a function of the traced Polyakov loop for G V = 0.25G s and for the canonical relation G V = 0.5G s found through Fierz transformation of the current-current interaction [19,20]. In order to construct these curves, we use Eq. (21) with the constituent quark masses found by the transcendental equations given by Eqs. (23) and (26). Φ is free to run in this case. Our analysis points out that in the range of chemical potentials from µ = M (vac) u,d to µ = Λ, it is possible to find curves in which Ω PNJL0 present global minima for nonvanishing values of the traced Polyakov loop, see Figs. 1c and 1d, for instance. For the curve presenting a 3 = −0.2 and G V = 0.25G s , we clearly see a single minimum of Ω PNJL0 at Φ = 0 for µ = M (vac) u,d (panel a). However, Φ ∼ 0.82 produces a minimum of Ω PNJL0 for µ = Λ (panel d). As in the SU(2) case [26], we verify a confinement/deconfinement phase transition in the model, i.e, the existence of null and nonvanishing values for the order parameter Φ as a function of the chemical potential. One can observe in panels e, f, g and h that a similar behavior is verified for the curves in which G V = 0.5G s . Another feature observed here, and also present in the SU(2) version of the model [26], is the need of G V = 0 values in order to become possible the existence of a global minimum of Ω PNJL0 at Φ = 0. For G V = 0, the only minimum is found at Φ = 0 for different a 3 values (figure not shown). The vector interaction still plays an important role even in the SU(3) version of the PNJL0 model.
The analysis of the Ω PNJL0 × Φ curves is also useful in order to determine the chemical potential in which the confinement/deconfinement phase transition takes place, named here as µ conf . For the model with G V = 0.25G s and a 3 = −0.1, for example, we see from Fig. 2 that µ conf = 533.15 MeV. This particular chemical potential value is defined as the one that generates a grand canonical potential presenting two minima with the same value for Ω PNJL0 , i. e., Ω PNJL0 (Φ 1 ) = Ω PNJL0 (Φ 2 ). In that case, Φ 1 and Φ 2 delimit the boundaries of the thermodynamical phases associated to quark confinement and deconfinement. The traced Polyakov loop is the order param-  eter of this transition. The procedure described above can be used to construct the Φ × µ curve depicted in gion determined by the condition given by µ < µ conf . On the other hand, deconfinement is achieved through a first order phase transition, in the region where µ > µ conf , according to the thermodynamics presented by the model. This feature was also observed in Ref. [37] in which the authors included the traced Polyakov loop in a hadronic SU(3) nonlinear σ model along with quark degrees of freedom. Furthermore, in Ref. [55] another version of a three-flavor PNJL model at zero temperature was studied. Specifically, a dependence on the quark density was introduced in the b 2 (T ) function of the Polyakov potential given in Eq. (11). In that model, a first order phase transition is also verified in the Φ(µ) curve. Another effect verified in Fig. 3 is the decreasing of µ conf as G V increases, exactly as in the SU(2) case [26]. Now we evaluate three-flavor PNJL0 model equations by solving Eqs. (23) and (26) together with condition ∂Ω PNJL0 /∂Φ = 0. In the SU(3) case, this procedure consists in searching for solutions of M u = M d , M s and Φ, obtained from the set of 3 coupled equations aforementioned. By following this method, we are able to evaluate Eq. (21), with result displayed in Fig. 4. One can notice a typical signal of a system presenting first order phase transition [56][57][58], since Fig. 4 shows that Ω PNJL0 (thermodynamical potential that describes the system) is multi-valued. A requirement from thermodynamical stability [56] imposes that branches associated to unstable and metastable solutions is removed from the final curve. The "crossing points" indicate the chemical potential values related to the first order phase transitions. In the case of the PNJL0 model, there are two of them, one from confinement/deconfinement and another one from restored/broken chiral symmetry of the light quarks, that give rise to µ conf = 533.15 MeV and µ chiral = 372.10 MeV, respectively, for the parametriza- Notice that at µ = µ conf , we find the same value for Ω PNJL0 related to the two global minima of Fig. 2, namely, Ω PNJL0 (µ conf ) ∼ −331 MeV/fm 3 . This feature confirms the equivalence of the two methods used to define the first order phase transition points in such systems. It is also worth to notice the decrease of µ conf as G V increases, exactly as pointed out in Fig. 3.
In the region where µ chiral < µ < µ conf we can recognize, as in the SU(2) version of the model [26], a particular thermodynamical phase in which the light quarks present a very low mass in comparison to their vacuum values (almost fully restored chiral symmetry), and with Φ = 0 (confined quarks), i.e, the so called quarkyonic phase, pointed out in Refs. [46,[59][60][61][62][63], for example. The quarkyonic phase related to the light quark sector can also be identified by looking at the constituent quark masses as a function of µ, as displayed in notice two sharp phase transitions occurring at µ chiral = 372.10 MeV and µ conf = 533.15 MeV, for G V = 0.25G s , with the difference defining the "size" of the phase, in this case, ∆µ = µ conf − µ chiral = 161.05 MeV. As we can also see, this difference is decreased for G V = 0.5G s . More specifically, for the SU(3) PNJL0 model ∆µ depends on both, the strength of the vector channel, G V , and the a 3 parameter related to the Polyakov potential. In Tables I  and II we present some values of ∆µ by running a 3 and G V independently. From Fig. 5, we can also notice a strong reduction of the constituent strange quark mass at µ = µ conf . At this point, M s undergoes a reduction of almost 50% for G V = 0.25G s , while for G V = 0.5G s the reduction is somewhat larger. The origin of this feature is the emergence of Φ exactly at this chemical potential, see Fig. 3. The decreasing of the couplings G s (Φ), G V (Φ) and K(Φ) when Φ becomes nonvanishing induces a reduction in the constituent quark masses, including the strange quark one. We reinforce that such an effect is not present in the original SU(3) NJL model, also show in Fig. 5 (in that case changes in M s are only gradual, see the dotted line). For the light sector, on the other hand, it is also verified a decreasing of M u,d in the PNJL0 model, however, with much smaller intensity, since at µ = µ conf the light quark masses are nearly vanishing already.
As already pointed out, the effect of the traced Polyakov loop depends on G V and a 3 . In Fig. 6 for a 3 = −0.27 and G V = 0.25G s , we verify that µ conf is equal to the cutoff parameter Λ resulting a smaller reduction of M s with respect to the NJL model, differently of what is seen in Fig. 5 in the previous case. For the sake of comparison, we also show the results concerning the parametrization presenting G V = 0.5G s (value obtained from the Fierz transformation of the current-current interaction). Notice that the decrease of the reduction in M s is also observed for this case.

C. Charge neutral system of quarks and leptons in weak equilibrium within hybrid stars
In this section we proceed to show how the confinement effects described by the traced Polyakov loop affect the charge neutral system of quarks and leptons (muons and massless electrons) in weak equilibrium. This study is relevant in the context of quark stars or hybrid stars [49,[64][65][66][67][68][69][70], constructed from a particular quark model at zero temperature and at high density regime. This is a system fulfilling the following conditions for the quarks and leptons chemical potentials: µ s = µ d = µ u + µ e , and µ e = µ µ . It is possible to write µ u , µ d and µ s in terms of a common quark chemical potential (µ), and in terms of the electron one (µ e ) as in which the relation µ B = 3µ holds. Furthermore, the requirement of charge neutral quark matter leads to the additional relation to be satisfied in which ρ e is the electron density. µ e and ρ e are related to each other through ρ e = µ 3 e /(3π 2 ). The muon density is ρ µ = [(µ 2 µ − m 2 µ ) 3/2 ]/(3π 2 ), and m µ = 105.7 MeV. In this case, the total energy density and total pressure are given, respectively, by and We start by showing in Fig. 7 the µ B dependence of the quark fractions defined as Y i = ρ i /(ρ u + ρ d + ρ s ), for the SU(3) PNJL0 model with G V = 0.25G s and different values of a 3 . We clearly verify that d quark fraction gradually reduces to a value close to the u quark fraction. Furthermore, at the chemical potential related to the deconfinement transition, point in which Φ becomes nonvanishing, one observes that this reduction is anticipated by a sharp variation. For the s quark fraction, the opposite is found, namely, Y s approaches to Y u with sharp variations produced by the emergence of Φ. In the case of Y u , we find no significant variations along the µ B dependence in comparison with both previous cases. From the figure, it is also clear that the critical baryonic chemical potential moves to the direction of increasing µ B as |a 3 | increases, feature also present in symmetric quark matter, as shown in Table I.
As mentioned before, an important application of a quark model available at zero temperature is the description of hybrid stars, in which a hadron-quark phase transition takes place at high density regime. A recent study based on a model-independent analysis, performed in Ref. [71], suggests evidences for massive stars composed by cores of quark matter. Here we take the PNJL0 model to take into account the quark side of the equations of state used to describe stellar matter. For the hadronic side, we use the relativistic mean-field model DDHδ given by the following Lagrangian density where for i = σ, ω, and (38) for i = ρ, δ. In Eq. (36) ψ is the nucleon field whereas σ, ω µ , ρ µ , and δ are the scalar, vector, isovector-vector, and isovector-scalar fields, that represents mesons σ, ω, ρ, and δ, respectively. The antisymmetric tensors F µν and B µν are written as F µν = ∂ ν ω µ − ∂ µ ω ν and B µν = ∂ ν ρ µ − ∂ µ ρ ν . The nucleon rest mass is M nuc , and the mesons masses are m σ , m ω , m ρ , and m δ . A more simplified version of this model was firstly proposed by J. D. Walecka in Ref. [72] in which the couplings between mesons and nucleons are constants. Here we choose to use a more sophisticated structure for the hadronic model in which such couplings are density dependent, namely, Γ i (ρ), with ρ being the sum of the densities of proton (ρ p ) and neutrons (ρ n ). In Eqs. (37) and (38) ρ 0 is the nuclear matter saturation density and Γ i (ρ 0 ), a i , b i , c i , and d i are free parameters of the model, that along with the nucleon and meson masses, define a specific hadronic parametrization. In particular, we use here the DDHδ [73] one. It is one of the 263 parametrizations tested against a set of constraints related to symmetric nuclear matter, pure neutron matter, symmetry energy and its slope (evaluated at ρ = ρ 0 ) [74]. It was also shown that this particular model is capable to generate massive neutron stars with mass around two solar masses [75].
Concerning the inclusion of more baryons in the hadronic sector, we remind the reader that some studies point out to a possible suppression of the hyperons population in hadron-quark phase transitions [69,76] due to the onset of hyperons occurring above the onset of quark matter. Furthermore, the inclusion of hyperons soften the hadronic equation of state with a direct consequence of producing not so massive stars in comparison with the case in which only nucleons are considered. Because of that, and since the main focus of this work is to analyze the thermodynamical structure of a quark model presenting deconfinement effects at zero temperature regime (PNJL0 model), we decide to use a RMF model only with nucleons, avoiding a proliferation of more coupling constants to be fixed in the hadronic sector (hyperons interactions are not completely determined and some assumptions are required).
The Euler-Lagrange equations are used in order to obtain the field equations of the model. In addition, the implementation of the Hartree approximation (no Fock term) [72,77] in these equations, leading to σ → σ ≡ σ, ω µ → ω µ ≡ ω 0 , ρ µ → ρ µ ≡ρ 0(3) , and δ → < δ >≡ δ (3) , allows the determination of the energy-momentum tensor, quantity used to calculate the energy density and pressure of the hadronic model. These last two thermodynamical equations of state are given by and ρ 3 = ρ p − ρ n , and ρ p,n = k F 3 p,n /3π 2 . Furthermore, the fields are determined as σ = Γ σ (ρ)ρ s /m 2 σ , ω 0 = Γ ω (ρ)ρ/m 2 ω ,ρ 0(3) = Γ ρ (ρ)ρ 3 /2m 2 ρ , and δ (3) = Γ δ (ρ)ρ s3 /m 2 δ with ρ s3 = ρ sp − ρ sn . The nucleon effective masses read with (−) for protons and (+) for neutrons. Finally, the proton and neutron chemical potentials, obtained through µ p,n = ∂E HAD /∂ρ p,n , are given by and respectively. We also implement charge neutrality and chemical equilibrium in the hadronic side by including muons and massless electrons in the system, in the same way it was done in the quark sector. This leads to the following relations: ρ p − ρ e = ρ µ and µ n − µ p = µ e , with µ µ = µ e . In this case, total energy density and total pressure are constructed as E tH = E HAD + E e + E µ and P tH = P HAD + P e + P µ . Moreover, E e,µ and P e,µ can be extracted from Eqs. (34) and (35). In Fig. 8a we display total pressure as a function of µ B for the hadronic model and for three different parametrizations of the PNJL0 quark model. In the hadronic side, one has µ B = µ n . For the quark sector we construct the PNJL0 parametrizations as follows: first we choose some values for the ratio G V /G s , namely, 0.15, 0.25, 0.35 and 0.5. For each of these values we impose that the point related to the confinement/deconfinement phase transition takes place exactly at the match with the hadronic equation of state. The value of a 3 is then determined for this purpose. This procedure, also adopted in Ref. [55], generates the linked models presented in Fig. 8b, in which the entire curves start with the hadronic model until the crossing point, and change to the deconfined quark model thereafter. This is the basis of the Maxwell construction, that imposes equal pressures and chemical potentials at both phases, for a fixed temperature, with a sharp first-order phase transition. The Notice that in set 4, where G V /G s = 0.5, we found a value very close (a 3 = −0.458) to one used in Ref. [37] (a 3 = −0.4). In Fig. 9 it is depicted how P t depends on E t for the linked models. It is clear that the effect of the first-order phase transition is to establish a plateau in the pressure with a gap in the energy density. The values of the transition pressure (P trans ) and the energy density gap at the transition point (∆E trans ) for each PNJL0 parametrization are also shown. Another interesting feature present in this analysis is that the two quantities, P trans and ∆E trans , increase with G V , or equivalently, with |a 3 |, both free parameters related to the PNJL0 quark model. The increasing of P trans means that the hadronic side of the linked model persists at higher values of µ B . This leads to a decreasing of the quark side, consequently, whilst the increasing of ∆E trans indicates an increasing of the coexistence region of hadrons and (deconfined) quarks. In order to describe a spherically symmetric hybrid star of mass M , we solve the Tolman-Oppenheimer-Volkoff (TOV) equations [78][79][80] given by (G = c = 1) whose solution is constrained to p(0) = p c (central pressure) and m(0) = 0. At the star surface, one has p(R) = 0 and m(R) ≡ M , with R defining its radius.
The equations of state used as input to solve Eqs.(46)-(47) are given by total pressure and total energy density of the linked models presented before. Furthermore, at the hadronic side we describe the crust of the star by the model developed by Baym, Pethick and Sutherland (BPS) [81] in a density region from ρ = 0.16×10 −10 fm −3 to ρ = 0.89 × 10 −2 fm −3 . The mass-radius profiles of the hybrid stars obtained from the models studied here are shown in Fig. 10. In this figure we compare our re- , respectively, and the recent one related to the MSP J0740+6620 pulsar at 68.3% credible level, namely, M = 2.14 +0.10 −0.09 M ⊙ [84]. As one can see, the hybrid model is able to produce stars with masses inside these boundaries. For the sake of comparison, we also depict in Fig. 10 the recent data extracted from the Neutron Star Interior Composition Explorer (NICER) mission concerning the radius related to the stars of masses equal to 2.08M ⊙ and 1.4M ⊙ . Such points are given by M = 1.4M ⊙ with R = (12.45 ± 0.65) km [85], M = (2.08 ± 0.07)M ⊙ with R = (12.35 ± 0.75) km [85]  . We also observe that hybrid star mass-radius curve is essentially dominated by the hadronic phase, leaving little room to the detailed dynamics from the PNJL0 model, unless for stars with the smallest radius, where the mass is somewhat decreasing by shrinking its size. That effect comes from the softening of the repulsion in the PNJL0 model when the deconfinement phase transition is approached, without the star loosing the stability, as we discuss below.
As an illustration, we plot in Fig. 11 how the pressure depends on the radial coordinate for a particular hybrid star constructed from the four sets of the DDHδ-PNJL0 model. In this case we chose the M = 2.098M ⊙ star. Notice that this particular star is predicted by all the sets considered here. In this figure, the kinks identify the quark cores of the star with radius R q ∼ 1.0 km for set 1. For sets 2 and 3 we found practically the same result of R q ∼ 1.5 km and R q ∼ 1.6 km, respectively. Finally, set 4 predicts a quark core of R q ∼ 1.3 km. In our approach, the core is composed by deconfined quarks.
On the other hand, hadrons are the degrees of freedom of the star from r = R q to r = R. As a last remark of our results concerning the massradius diagrams, we observe that for each parametrization constructed, namely, sets 1 to 3 of the DDHδ-PNJL0 hadron-quark model, one identifies linear branches in which M decreases as R decreases. From the investigation of the static condition for the star's stability, this kind of branch would indicate unstable configurations [88,89]. However, a dynamical analysis of the star stability can give further insights by verifying its response to small radial perturbations, see for instance Refs. [91][92][93][94][95][96][97] for details. By using this criterion, stars are said stable if such perturbations produce well defined oscillations. On the other hand, an indefinite increasing of the perturbations amplitude characterizes unstable stars. The set of differential equation that need to be solved, coupled to the TOV ones, are given by [91][92][93][94][95][96][97] and d∆p dr = ξ ω 2 e λ−ν (p + ǫ)r − 4 dp dr + ξ dp dr 2 r p + ǫ − 8πe λ (p + ǫ)pr + ∆p dp dr with e λ = 1 − 2m(r)/r, dν/dr = −2(dp/dr)(p + ǫ) −1 , and Γ = (ρ/p)(dp/dρ). The relative radial displacement and the perturbation of the pressure, ξ = ∆r/r and ∆p, respectively, are assumed to have a time dependence of e iωt with ω being the eigenfrequency. Furthermore, other important aspect that must be mentioned is the kind of phase transition chosen. In stars with two phases, transitions are classified as "fast" or "slow" depending on the timescale of the transition compared to the timescale of the fundamental oscillation [90]. As we can see in Refs. [92][93][94], slow transitions have important consequences in the microscopic structure of compact stars with two different phases. It is well known that in stable stars the fundamental frequency has to satisfy the condition of ω 2 > 0, and the last stable star occurs at ω = 0. In general, this property coincides with the maximum mass in the mass-radius diagram. However, when a star presents two phases and the transition is slow, the last stable star (ω = 0) can occur after the point of maximum mass, as we can verify in the results shown in Ref [92], for instance. By following this approach, we verify that all branches analyzed in the parametrizations shown in Fig. 10 are stable, i.e., we found no points in which ω = 0.

IV. SUMMARY AND CONCLUDING REMARKS
In this work we generalize the study performed in Ref. [26], where we have proposed a version of the Polyakov-Nambu-Jona-Lasino that exhibits effects of the confinement/deconfinement phase transition at zero temperature regime for the two-flavor case. We improve our analysis in order to take into account the strangeness in the dense system by introducing the dynamics of the strange quark s, whose current mass is considered as m s = 140.7 MeV. We start with the Nambu-Jona-Lasinio Lagrangian density model with the 't Hooft interaction included. We recall the NJL grand canonical potential written in terms of the strengths of the scalar, vector and quark mixture channels, namely, G s , G V and K, respectively. The implementation of the infrared quarkgluon dynamics is performed by making these couplings dependent on the traced Polyakov loop as follows, G s → G s (G s , Φ) = G s (1−Φ 2 ), G V → G V (G V , Φ) = G V (1−Φ 2 ), and K → K(K, Φ) = K(1 − Φ 2 ). The motivation for this procedure is to ensure vanishing quark interactions at Φ = 1 (deconfined phase). By using these new functions G s (K, Φ), G V (K, Φ) and K(K, Φ), it is possible to define a new Polyakov potential now with strangeness included in the back reaction (quarks and gluons affecting each other).
We also analyze how the confinement/deconfinement phase transition takes place from the study of the minima of Φ at a fixed (common) quark chemical potential. In this investigation, we focus on the case of symmetric quark matter, namely, the system presenting µ u = µ d = µ s . The same study applied to other values of µ generates a typical first order phase transition for the order parameter Φ as a function of µ. Therefore, the system presents two phase transitions, namely, the one defined by Φ, and the other related to the broken/restored chiral symmetry phase transition, in which the quark condensates are the order parameters. This last one is already presented in the original NJL model. The region between these two transitions can be identified, as in the SU(2) case, as the quarkyonic phase, defined here as the phase in which the light quarks present restored chiral symmetry but are still confined (Φ = 0).
Another interesting feature of the SU(3) version of the PNJL0 model is the strong reduction of the constituent strange quark mass (M s ) at the deconfinement phase transition. We verify that M s undergoes a reduction of almost 50% for the parametrization in which G V = 0.25G s and a 3 = −0.1. The dynamics of the traced Polyakov loop at T = 0 favors the restoration of chiral symmetry even for the strange quark. This effect is not present in the original SU(3) NJL model since Φ is always vanishing in this case.
We also implement the chemical equilibrium and charge neutrality in the SU(3) PNJL0 model in order to apply it in the description of hybrid stars. For the hadronic side we use a density dependent model and match both equations of state exactly at the point where quarks are deconfined. Different PNJL0 parametrizations were constructed by changing the values of G V and a 3 . We verify that the increasing of both, G V and |a 3 |, increases the transition pressure and the gap in the energy density of the hadron-quark phase transition. Finally, we generate the mass-radius profiles determined by the parametrizations of the PNJL0 model (different G V and |a 3 | values) coupled to the hadronic model. Our results indicate that it is possible to construct hybrid stars compatible with some observational data. In particular, the recent ones proposed by the NICER mission [85][86][87].