1 Introduction↩︎

The theoretical description of strongly interacting matter at high baryonic densities plays a central role in investigations aiming to describe the structure of neutron stars (NS). In general, hadrons and/or quarks are considered to be the most relevant degrees of freedom representing the matter that makes up these dense stellar objects. In particular, the possibility that hadron and quark matter form a hybrid NS became more appealing due to the recent observational data of two solar mass pulsars [1]–[3] which suggest that quark matter may indeed be present within their core [4], [5]. In this case, an appropriate equation of state (EoS), preferably obtained from the fundamental quantum chromodynamics (QCD) theory, should be used to describe the ultra-dense core region.

An ab-initio determination of the QCD EoS for cold and compressed quark matter represents a challenging task, since lattice QCD (LQCD) calculations are limited to low baryon chemical potential values (\(\mu_B\)) due to the well-documented sign problem [6], which obstructs extensions to intermediate and large \(\mu_B\) values. In the low baryon-density regime, chiral perturbation theory (chEFT) provides a precise description [7], [8]. At intermediate and higher \(\mu_B\) values, one can consider effective models [9]–[11] as well as perturbative QCD (pQCD) [12]–[16], but unfortunately these alternatives also present some shortcomings. For instance, most effective quark models (such as the Nambu-Jona-Lasinio model [17]–[20] and the quark-meson model [21]) do not take asymptotic freedom and gluonic degrees of freedom into account. At the same time, pQCD applications are reliable only at high baryon densities of the order of \(\rho_B \sim 25-40 \, \rho_0\) (\(\rho_0 = 0.16 \, {\rm fm}^{-3}\)) [22], where asymptotic freedom allows for weak coupling expansions. However, for NS with masses \(\sim 1.4 − 2M_\odot\) the typical central densities are expected to be around the \(5 − 10\, \rho_0\) range [23] where quark matter is still strongly coupled. A bridge between the low- and high-density regions was elaborated via model-agnostic equations of state (EoS) of compact stars applying different statistical methods (see, e.g., [4], [24]–[29]), via thermodynamical constrains[30], and via microscopic models within a Bayesian inference approach (see, e.g., [31]–[37]), imposing NS observational constraints in addition to those from ab-initio pQCD and chEFT calculations. Such approaches are remarkable in that they produce EoSs over the largest possible \(\mu_B\) range. However, they are limited by the present accuracy of theoretical chEFT and pQCD predictions at both ends and by the accuracy of available data. This is also complicated by the phase transition [38] at intermediate \(\mu_B\) values, which generally requires further model-dependent assumptions.

Concerning higher order pQCD evaluations, it is important to recall that these are notoriously complicated by thermal or medium properties [39]–[41]. In particular, infrared divergences proliferate starting at order \(\mathcal{O}(\alpha_s^2)\), which cancel only through a specific resummation of a specific class of Feynman diagrams, but producing non-analytical contributions to the naive weak-coupling perturbative expansion. For cold quark matter, these contributions are classified as “soft”, associated to the scale \(m_E\sim \sqrt{\alpha_s}\mu\), in contrast to the “hard” contributions provided by quarks at the Fermi surface given by the quark chemical potential scale, \(\mu\). The groundbreaking calculation of Freedman and McLerran [12] revealed the emergence of an \(\alpha_s^2\ln\alpha_s\) dependence (in the massless quark approximation) in the next-to-next-to-leading order (NNLO) pressure. These soft logarithm contributions emanate from the plasmon (ring) resummation of soft infrared divergences, nowadays well-understood also within the hard thermal loop (HTL) [42] effective field theory framework [16], [43], which has also led to the all-order resummation of the leading soft logarithms [44]. The cold quark matter pressure calculations have been generalized to include quark masses [13], [45]–[47]; and thermal effects [48], [49]. In the recent few years, even the full NNNLO is in sight [16], [22], [43], [50] in the massless quark approximation, expecting significantly improved pQCD accuracy. Despite this remarkable recent progress, at lower \(\mu_B\) beyond the perturbative regime, the pQCD pressure exhibits a loss of accuracy, mainly reflected in a rapidly growing sensitivity on the arbitrary renormalization scale \(\Lambda\) in the (\({\overline{\rm MS}}\)) scheme[13]. Although the pQCD NNLO pressure is perturbatively renormalization group (RG) invariant, this implies that while formally of higher order \({\cal O}(\alpha_s^3)\), the residual scale dependence has sizable impact, even for not that small \(\mu_B \sim 2 \, {\rm GeV}\), which corresponds to a baryon density \(\sim 40\, \rho_0\). It indicates that possible non-perturbative effects remain quite important in the lower/intermediate density range relevant to NSs. Moreover, it appears that this large residual scale dependence originates principally [43], [44] from the hard sector, contributing predominantly to the total pressure, while in comparison the soft sector is subdominant and under better control [22], [43], [44]. This is further worsened if accounting for the strange quark mass effects, not always negligible at moderate and low \(\mu_B\) values. For massive quarks, pQCD calculations are even more involved and the state-of-the-art is limited to NNLO [13], thus impinging on scale uncertainty reductions.
In the present work, we consider a resummation approach which incorporates higher order RG properties in order to mitigate the above mentioned pQCD problems in evaluations aiming to describe the quark sector of NSs. The renormalization group optimized perturbation theory (RGOPT) [51]–[53] is based on a modified perturbative expansion, around massive quasiparticle states regularizing infrared divergences. Basically, this is similar to other approaches at finite temperature and densities based on interpolated Lagrangians augmented by prescriptions aiming to go beyond strictly perturbative expansions. Such screened perturbation theory [54]–[56] (SPT) for the thermal \(\phi^4\) model has found many applications, especially when used within the HTL Lagrangian, hence known as Hard Thermal Loop perturbation theory (HTLpt) [57], [58]. A specific feature of the RGOPT prescription is that it generates a “RG-dressed” screening mass, entailing an all-order RG-driven \(\alpha_s\) dependence. At finite temperature, RGOPT has been applied to \(\phi^{4}\) theory [59] up to NNLO [44], and to NLO for the hot QCD pressure [60], [61], where in both cases it drastically reduces the residual renormalization scale dependence with respect to the standard weak-coupling expansion, SPT or HTLpt. For cold quark matter, the NLO RGOPT [62] also reduces the residual scale dependence, although more moderately than in the \(T\ne 0\) case. Recently, it has been also used to describe non-strange and strange quark stars (QS) at NLO [63], [64], yielding results which are less sensitive to scale variations than those provided by pQCD.

Meanwhile, the method has been considered at NNLO to describe the \(N_f=2+1\) QCD EoS in the case of symmetric matter [65], \(\mu_u=\mu_s=\mu_d\equiv \mu\). In the present work, the NNLO results of Ref [65] will be extended to treat non-symmetric quark matter in \(\beta\)-equilibrium, so that an EoS tailored to describe pure quarks stars as well as the core of hybrid NSs can be made available. For practical purposes, we will also provide an approximate but accurate compact formula for our obtained (thermodynamically consistent) pressure as a function of the renormalization scale and baryon chemical potential, so that the EoS can be readily computed. When applying pQCD to cold and dense matter, one usually sets \(\Lambda = X \, \mu_B/3\) where, for \(N_f = 2+1\) flavors, \(\mu_B\) relates to the individual quark chemical potentials through \(\mu_B = (\mu_u+\mu_s+\mu_d)\). In general, \(X\) is probed for \(1 \le X \le 4\) with \(X=2\) representing a “central" fiducial scale.

In the present study, we will compare state-of-the-art NNLO pQCD and RGOPT predictions including massive quarks, also aiming to constrain the renormalization scale, by imposing agreement with NS observational data as in Refs. [63], [64]. In particular, we will consider the mass and radius measurements from NASA’s Neutron Star Interior Composition Explorer (NICER) mission for the pulsars PSR J0030+0451 [66], [67], PSR J0740+6620 [68]–[71], PSR J0437+4715 [72], PSR J0614-3329 [73], and the light compact object HESS J1731-347 [74]. A quite strong constraint is imposed by the pulsar J0952−0607, a “black widow” pulsar with a gravitational mass above 2.27\(M_\odot\) (2.12\(M_\odot\)) at 1\(\sigma\) (3\(\sigma\)) confidence interval [75], [76]. However, the mass measurements of pulsars as PSR J0952−0607 depend on the optical modeling, which introduces larger uncertainties than estimates from pulsar timing, as done by NICER observations. Additional constraints are imposed by the detection of gravitational waves by the LIGO Virgo Collaboration, in particular GW170817 [77], [78] from a binary NS merger and GW190814 from the coalescence of a black hole with a 2.5-2.67 solar mass compact object, possibly a NS. Apart from discussing pure QSs, we will also build some hybrid EoSs by considering a Maxwell construction to describe a first order phase transition between the hadronic and quark phases, obtained within the NNLO RGOPT. In order to cover the present existing uncertainty in the hadronic phase, we select a few nucleonic models representative of soft and stiff equations of state to build the hybrid EoS, all based on a relativistic mean-field description of nuclear matter: SFHo [79] and DIDY [80] which represent a soft EoS and DD2, which can produce a reasonably hard EoS [81]. These equations of state satisfy the well-known nuclear matter and NS properties.

The work is organized as follows. In Section 2, the quark pressure \(P(X,\mu_B)\) is determined within the RGOPT resummation framework (this section may be skipped by readers mainly interested in applications to compact stars). In Sec. 3 we present pocket formulas for the thermodynamically consistent pressure, for the case of \(\beta\)-equilibrium and charge neutral matter. The RGOPT and pQCD results for quark and hybrid stars are presented and compared in Sec. 4. Finally, in Sec. 5 we present our conclusions.

2 RGOPT pressure↩︎

In this Section, we review the derivation of the RGOPT pressure for cold dense quark matter at NNLO by summarizing the derivation contained in Ref. [65]. As mentioned in the Introduction, by construction the RGOPT aims at reducing the residual renormalization scale dependence as compared to the standard weak coupling expansion. Concerning specifically the cold and dense pQCD pressure, originally derived at NNLO with full quark mass dependence [13], as emphasized above, its sizable residual renormalization scale dependence (thus originating from unknown higher perturbative orders) originates predominantly from the hard sector contributions, which is further enhanced for massive quarks [65]. Accordingly, we consider the RGOPT only in the massive quark sector for simplicity¹.
Since this approach is basically a modification of the standard weak coupling expansion of a massive theory, we start from the NNLO weak coupling expansion in \(\alpha_s=g_s^2/(4\pi)\) of the cold and dense quark matter pressure, whose full NNLO quark mass dependence was originally calculated in Ref. [13]. The relevant Feynman graphs up to NNLO in the weak coupling expansion are given in Figs. 1 and 2.

Importantly, instead of considering massless \(u,d\) and one massive strange quark as usual in pQCD, we anticipate that the RGOPT prescription induces an (all-order) dressed screening quark mass \(m ={\cal O}(g_s\, \mu)\) common to the three flavors, as explained below (while we still approximate vanishing current \(u,d\) quark masses). Thus, all graphs in Figs. 1, 2 are evaluated with massive quarks. Moreover, non-diagonal contributions occur from ultimately distinct quark masses \((m+m_u,m+m_d,m +m_s)\), where \(m_s\) is the strange quark current mass and \(m_u=m_d\equiv 0\). The graphs in Fig. 2 produce such non-diagonal terms at NNLO, derived in Ref. [65], to which we refer for more details².

The resulting NNLO pressure for generic quark masses \(m +m_i\) (that we refer to in the sequel as \(P^{N_f=2^*+1^*}\) to distinguish from standard \(N_f=2+1\) pQCD pressure), with quark chemical potential \(\mu_i\), may be formally written as \[\begin{align} \label{Ptotal} &&P(m,m_i,\mu_i)= \nonumber \\ && \sum_{i=1}^{N_f}\left( P_{\rm LO}^{m}(m+m_i,\mu_i) +P_{\rm NLO}^{m}(m+m_i,\mu_i) +P_{\rm 2GI}^{N_f=2^*+1^*}(m+m_i,\mu_i) +P_{\rm VM}^{N_f=2^*+1^*}(m+m_i,\mu_i)\right) +P_{\rm Ring}^{N_f=3^*}\!(m) \nonumber \\ && +\sum_{i=1}^{N_f}\left(P_{\rm LO}^{v}(m+m_i) +P_{\rm NLO}^{v}(m+m_i) +\sum_{i=1}^{N_f-1} \left(P^{v}_{\rm NNLO,1}(m)+P_{\rm sub,1}(m)\right) +P^{v}_{\rm NNLO,2}(m+m_s) +P^{v,nd}_{\rm NNLO}(m,m+m_s)\right)\nonumber \\ && +P_{\rm sub,2}(m+m_s) +P_{\rm sub}^{nd}(m,m+m_s)\,, \end{align}\tag{1}\] where the different expressions may be found in Ref. [65] (see also Appendix 6 for a summary) and we briefly comment on these here. Note that up to NLO, the contributions factorize per flavor, while at NNLO non-diagonal terms occur, as indicated by the “nd” quantities in Eq. (1 ), due to the occurrence of independent quark loops in Fig. 2.

Accordingly, the NNLO contributions are conveniently decomposed into four classes:

1. The NNLO weak coupling expansion of the matter (\(m\)) contributions, \(P^{Nf=2^*+1^*}_{\rm 2GI}(m_i,\mu_i)\) in Eq. (1 ) corresponds to all two-gluon irreducible contributions in Fig. 1.

2. \(P_{\rm VM}^{N_f=2^*+1^*}(m_i,\mu_i)\) corresponds to mixed vacuum-matter (VM) contributions as indicated in Fig. 2.

3. The plasmon sum “ring” contributions, that result after the all-order resummation sketched in the rightmost graph of Fig. 2: \[P_{\rm Ring}(m) = -\Omega_{\rm Ring}^{N_f=3^\star}(m) \,.\] Note that the latter contribution could not be handled analytically while maintaining an exact dependence on arbitrary masses and chemical potentials. However, this \({\cal O}(\alpha_s^2)\) contribution is numerically an order of magnitude smaller than the NNLO \(P_{\rm 2GI+VM}^{N_f=2^*+1^*}(m,\mu)\). Therefore, we neglect any mixing altogether within this contribution, i.e. we take a universal quark mass \(m\) for the three flavors [65] (the actual \(m\) value being subsequently determined from the RGOPT prescription, as explained below). Contributions 1)–3) were originally derived in Ref. [13] with exact quark mass dependence, and extended to quark flavors with different masses in Ref. [65].

4. Finally, a fourth class of contributions, not present within the standard weak coupling expansion[13] and crucially involved in the RGOPT prescription, is to (re)introduce vacuum contributions to the pressure, \(P^{v}_{\rm LO}-P^{v}_{\rm NNLO}\) in Eq.(1 ), given by the \(\mu=0\) contributions in Fig. 1 and “VV” contribution in Fig. 2. Importantly, such vacuum contributions should be supplemented with zero-point “subtraction” contributions [53], [62] (represented in Eq. (1 ) by \(P_{\rm sub,i}\) terms) directly related to the vacuum energy anomalous dimension [82], [83], \(\propto m^4\,\mathbb{1}\), which mixes with the \(m \bar q q\) operator (see also Appendix 7). Often disregarded in the thermal literature on account of being medium-independent, those vacuum contributions play a crucial role for RG properties of a massive theory, as they are necessary to restore perturbative RG invariance of a massive pressure (equivalently vacuum energy)³.

Next, the actual RGOPT prescription is implemented by the following steps (see, e.g., Refs. [51], [62] for a review):

• First, the NNLO pressure in Eq. (1 ), being formally a function \(P(m,m_i,g_s^2,\mu_i)\), is modified by the following substitution (performed after standard coupling and quark mass renormalization): \[\{m_i, m \} \to \{ m_i, \eta\left(1-\delta\right)^a \}, \;\; g_s^2\to \delta g_s^2\,, \label{delta}\tag{2}\] where the exponent \(a\) is specified just below and, as already emphasized, the actual physical masses considered here are \(m_u=m_d=0, m_s \ne 0\). Initially, \(\eta\) is an arbitrary trial mass, thus distinct from the physical quark mass \(m_i\). Eq. (2 ) is such that \(0\le \delta \le 1\) interpolates between the massive (free) theory and the original theory respectively, with \(\delta\,\eta\) treated as an interaction. The pressure is then re-expanded in powers of \(\delta\) at the relevant perturbative order in \(g_s^2\) (here NNLO), and taking the limit \(\delta\to 1\) afterward to recover the original theory. Truncating the series at finite \(\delta\)-order in this way leaves a remnant \(\eta\)-dependence, fixed by RG properties as specified below. The latter prescription applied to Eq.(1 ) gives formally an expression denoted \(P^{N_f=2^*+1^*}_{\rm RGOPT}(\eta,g_s^2,m_i,\mu_i)\) below.

• A crucial feature of RGOPT is that the exponent \(a\) in Eq. (2 ) is uniquely fixed by the (massless) RG equation at LO to the critical value [51], [62] \[\label{acrit} a= \frac{\gamma_0}{2b_0}\;\; = \frac{4}{9} \Big|_{N_c=N_f=3} \,,\tag{3}\] where \(b_0\) and \(\gamma_0\) are the leading coefficients for, respectively, the beta function and the quark mass anomalous dimension (see Appendix B for notation). This prescription remains consistent with standard renormalization and maintains perturbative RG invariance⁴ even after the \(\delta\)-expansion is performed.

• Lastly, in related optimized perturbation (OPT) approaches (for instance, in SPT[55] and HTLpt at NLO[84], [85]) a standard prescription to fix the arbitrary mass parameter \(\eta\) is a stationarity principle, which in our present context takes the form: \[\label{eq:OPTequation} {\rm f}_{\rm OPT}(\eta, g_s,\cdots) \equiv \frac{\partial P^{N_f=2^*+1^*}_{\rm RGOPT}}{\partial \eta}\Big|_{\eta=\tilde{\eta}}=0\,.\tag{4}\] However, the number of solutions of Eq. 4 increases with perturbative order and the solutions are not necessarily real-valued due to nonlinear \(\eta\)-dependence. In contrast, the RGOPT prescription in Eq. 3 crucially guarantees that, at successive orders, a single solution \(\eta\) satisfies the leading behavior of the parameters RG flow in the original (massless) theory [51] (i.e. for QCD the solution satisfies asymptotic freedom). At the same time, to fix the trial mass parameter \(\eta\) beyond LO it appears more compelling to use the RG equation \[\label{eq:RGred} {\rm f}_{\rm RG}(\eta, g_s,\cdots) \equiv \left(\Lambda\frac{\partial}{\partial\Lambda}+\beta(g_s^2)\frac{\partial}{\partial g_s^2} \right) P^{N_f=2^*+1^*}_{\rm RGOPT}(\eta,g_s^2,\cdots)\Big|_{\eta=\tilde{\eta}}=0,\tag{5}\] with the beta-function \(\beta(g_s^2)\) (see Appendix B), instead of Eq. 4 (here for the simplified situation of vanishing physical mass \(m_i\)). Eq. (5 ) produces a “RG-dressed” screening quark mass, \(\eta(g_s^2,\mu)\), incorporating an all-order summation of certain topologies [62], [65].

Incidentally, for the cold quark matter pressure at NNLO, only Eq. 5 gives a solution which is perturbatively consistent with a screening behavior [65], \(\tilde{\eta}_{RG} = {\cal O}(g_s \mu)\) for \(g_s\to 0\), in contrast to Eq. 4 . Accordingly, our prescription is essentially uniquely defined at NNLO. Yet, such a RG solution is not guaranteed to be real-valued, as mentioned previously. Indeed, in the \({\overline{\rm MS}}\) scheme the above prescription for the quark matter pressure leads to a complex-valued \(\tilde{\eta}\) already at NLO [62]. We stress that this is a nonphysical artifact of the nonlinearity in \(\eta\) arising from insisting on an “exact” solution \(\tilde{\eta}(g_s,\mu)\). Alternatively, \(\tilde{\eta}(g_s,\mu)\) may be re-expanded perturbatively, which makes it always real-valued, but this generally loses [62], [86] part of the sought RG improvement from higher orders. To preserve the RG resummation benefits while circumventing the nonreal-solution issue, we rather exploit the renormalization scheme freedom to slightly shift perturbative coefficients from the original \({\overline{\rm MS}}\) ones, so as to restore a real-valued \(\tilde{\eta}\) solution. Our renormalization scheme change (RSC) is defined perturbatively [62] consistently prior to the modifications from Eq. (2 ), using \[\label{eq:RSC} \left(m, m_s\right) \to \left[m (1+B_2 g_s^4), m_s(1+B_2 g_s^4)\right]\,.\tag{6}\] After the modification induced from Eq. (2 ) and NNLO \(\delta\)-expansion, the arbitrary RSC parameter \(B_2\) is uniquely fixed upon requiring the minimal departure from \({\overline{\rm MS}}\): it can be shown that this amounts to solve for \(B_2\) and the trial mass \(\eta\) simultaneously the RG Eq. (5 ) and the constraint [51], [62] \[\label{eq:determinant} \rm det(f_{RG},f_{OPT})_{(g_s^2,\eta)} = \left(\frac{\partial f_{\rm RG}}{\partial g_s^2}\frac{\partial f_{\rm OPT}}{\partial \eta}- \frac{\partial f_{\rm OPT}}{\partial g_s^2} \frac{\partial f_{\rm RG}}{\partial \eta}\right)\Big|_{\eta=\tilde{\eta},\;B_2=\tilde{B}_2}=0 \,.\tag{7}\] Actually, when accounting for non-zero physical strange quark mass, we modify the RG operator, Eq. (5 ), as \[\label{eq:RGredms} \Lambda\frac{d}{d\Lambda}=\Lambda\frac{\partial}{\partial\Lambda}+\beta(g_s^2)\frac{\partial}{\partial g_s^2}-m_s\gamma_{m_s}(g_s^2) \frac{\partial}{\partialm_s} \,,\tag{8}\] with \(\gamma_m(g_s^2)\) the anomalous quark mass dimension (see Appendix B). Note that the genuine physical mass \(m_s\) thus remains unaffected by the \(\delta\)-expansion.
Applying the above RGOPT prescription to the cold dense \(N_f=2^*+1^*\) NNLO pressure Eq. (1 ) with respective flavor masses \((m,m,m+m_s)\), we obtain the doublet of real solutions (\(\tilde{\eta}\),\(\tilde{B}_2\)) and the RG-dressed mass spectrum \((\tilde{\eta}, \tilde{\eta}, \tilde{\eta}+m_s)\) respectively for \(u,d,s\) quarks. Accordingly, the resulting NNLO RGOPT pressure reads formally \[\label{Prgopt95nnlo} P_{\rm RGOPT} \equiv P^{N_f=2^*+1^*}_{\rm RGOPT}(\tilde{\eta}, g_s^2, \tilde{B}_2,\mu_i,m_s) \,,\tag{9}\] where we simply indicate its final dependence on relevant parameters since its explicit expression is quite involved [65]. With the aim of providing a simpler procedure for further applications, it is more convenient to construct a sufficiently accurate approximation through an appropriate fitting procedure, as described in the Section below.

3 Thermodynamically consistent RGOPT pressure in \(\beta\)-equilibrium↩︎

This Section details the construction of the \(\beta\)-equilibrated and thermodynamically consistent RGOPT pressure entering the EoS employed in the compact star applications of Section 4. Starting from the pressure in Eq. (9 ) (for symmetric matter) obtained upon applying the RGOPT prescription of Sec. 2, chemical equilibrium and charge neutrality are enforced by requiring the following relations to be satisfied \[\label{eq:BetaEq} \mu_u=\mu_d-\mu_e, \;\;\;\;\mu_s=\mu_d\,,\tag{10}\] \[\label{eq:charge95neutrality} {\rm f}_{\rm CN}\equiv\frac{2}{3}\rho_{\rm up}-\frac{1}{3}\rho_{\rm down}-\frac{1}{3}\rho_{\rm strange}-\rho_{\rm electron}=0\,,\tag{11}\] where \(\rho_i= \frac{\partial P_{\rm RGOPT}}{\partial\mu_i} |_{\Lambda}\) represents the number density for the \(i\)-th particle while \(\mu_e\) represents the electron chemical potential. In principle, thermodynamic quantities such as pressure and number density receive electronic contributions \(\sim \mu_e^4\), and \(\sim \mu_e^3\), respectively. Explicit checks confirmed that \(\mu_e \ll \mu\); therefore, these contributions are negligible, and we have subsequently omitted them from our analysis. Notice also that the determinant in Eq. (7 ) for the symmetric matter case is further generalized to incorporate the additional physical constraint of charge neutrality or, equivalently, positive \(\mu_e\). This is accounted for by a simple trick: substituting \(\mu_e\equiv\mu_A^2\) and enforcing the real-valuedness of \(\mu_A\) by extending Eq. (7 ) to: \[\label{eq:GeneralizedDeterminant} \rm det(\rm f_{RG},f_{OPT},f_{CN})_{(g_s^2,\eta,\mu_A)}=0\,,\tag{12}\] where “\(\rm f_{CN}\)” stands for the charge neutrality Eq. (11 ).

Note that, in the resulting pressure, all parameters (except \(\mu_i\)) exhibit a highly non-trivial dependence on the renormalization scale \(\Lambda\): manifestly so for \(\alpha_S(\Lambda),~ m_s(\Lambda)\), dictated by RG running, but equally for the additional parameters \(\tilde{\eta}, \tilde{B}_2\) that entail a non-trivial \(\Lambda\) dependence from the RGOPT prescription. Thermodynamic consistency is therefore spoiled, since the relation \(d P_{\rm RGOPT}/d\mu_i = \rho_i\), does not hold [13], [63], [87]–[89]. Nevertheless, one can obtain a thermodynamically consistent (TC) pressure by considering \[\label{Pth} P_{\rm RGOPT}^{TC}(...,\mu_i) \equiv P_{\rm RGOPT}(...,\mu_i) + \int_{\mu_0}^{\mu_i} d\mu \left( \frac{\partial}{\partial_\mu} P_{\rm RGOPT}(\mu) - \frac{d}{d\mu} P_{\rm RGOPT}(\mu) \right) \,,\tag{13}\] where \(\mu_0\) is specified below. Note that the additional term in Eq. (13 ) is formally of higher perturbative order: by construction, our NNLO \(P_{\rm RGOPT}\) is RG invariant up to neglected higher order, i.e. \(d P_{\rm RGOPT}/d\Lambda = {\cal O}(\alpha_s^3)\), then from \((d/d\mu) P_{\rm RGOPT} = \partial_\mu P_{\rm RGOPT} + (d\Lambda/d\mu) d P_{\rm RGOPT}/d\Lambda)\), the integrand above is \((d\Lambda/d\mu)\times {\cal O}(\alpha_s^3)\), and conventionally, one chooses a scale \(\Lambda = {\cal O}(\mu)\) within a certain range, so that \((d\Lambda/d\mu)\) is just a number. Adopting a thermodynamically consistent prescription for the compact star EoS is nonetheless essential, not only in general, but also within our construction, where Eq. (9 ) (and its generalization for non-symmetric matter) captures a non-trivial dependence on higher orders. The initial value \(\mu_0\) in Eq. (13 ) is chosen such that the total quark number density vanishes, \(\rho_{tot}(\mu_0)=0\).
Thus, by solving simultaneously Eqs. (12 ) and (8 ), we obtain in this case the triplet of solutions \((\tilde{\eta},\tilde{B}_2, \tilde{\mu}_A\)) (hence \(\tilde{\mu}_e^2\)) leading to a fully determined pressure in Eq. (9 ) for a given quark chemical potential \(\mu_i\) and renormalization scale \(\Lambda=X \mu_B/3\).

3.1 Compact formula↩︎

In order to reproduce more simply the quite involved NNLO RGOPT results, in this subsection we present fitting functions for the thermodynamically consistent NNLO RGOPT pressure given in Eq. (13 ), which are solely functions of the scale parameter \(X= 3\Lambda/\mu_B\) and the baryon chemical potential \(\mu_B\). Constructing a fitting function that accurately reproduces our results for all relevant \(\mu_B\) and \(X\) values is quite challenging. In practice, we conveniently provide two different formulas in two different \(X\)-regimes, \(3\le X\le 6\), which corresponds to the range essentially relevant for QS, and \(2\le X\le 3\), corresponding to the relevant range for hybrid stars. In general, values of \(X<2\) tend to match typical hadronic EoSs at very large \(\mu_B\). Consequently, the transition to quark matter occurs at chemical potential values that lie beyond the stability limit. Therefore, by not producing a quark core, such configurations can be safely neglected in the present application.

3.1.1 Current quark mass and QCD coupling input↩︎

In our numerical applications, for the running coupling \(g_s^2(\Lambda)\equiv 4\pi\alpha_s(\Lambda)\), we use the exact NLO result obtained for a given renormalization scale \(\Lambda\) upon solving \(g_s^2\) exactly from \[\label{eq:RunningNNLO} \Lambda_{\rm \overline{MS}}= \Lambda e^{-\frac{1}{2b_0 g_s^2}}\left(\frac{b_0 g_s^2}{1+\frac{b_1}{b_0} g_s^2} \right)^{-\frac{b_1}{2b_0^2}} \,,\;\tag{14}\] and fixing \(\Lambda_{\rm \overline{MS}}=330\;\)MeV [90] so that \(\alpha_s(\Lambda = 1.5 \, {\rm GeV})\simeq 0.326\) [91]. Since the expected RG scale dependence cancelations always occur between running coupling (or masses) at order \(\alpha_s^k\), and explicit \(\ln \Lambda\) dependence at order \(\alpha_s^{k+1}\), the NLO running is sufficient in principle when considering an NNLO pressure. Using higher order running coupling would hardly give any visible differences in our results below for the relevant range of parameters considered.

We recall that, in all applications below, we approximate the current masses \(m_u, m_d\) to zero, as usual. Concerning the strange quark, since mass renormalization is only needed at NLO \({\cal O}(\alpha_s^2)\) in our case, we use the NLO running mass, given in our normalization as: \[\label{msrun} m_s(\Lambda)=m_s(\Lambda_0)\left(\frac{g_s^2(\Lambda)}{g_s^2(\Lambda_0)}\right)^{\frac{\gamma_0}{2b_0}}\left(\frac{1+\frac{b_1}{b_0}\,g_s^2(\Lambda)}{1+\frac{b_1}{b_0}\,g_s^2 (\Lambda_0)}\right)^{\frac{\gamma_1}{2b_1}-\frac{\gamma_0}{2b_0}} \,,\tag{15}\] where \(m_s(\Lambda_0=\mathrm{2 GeV}) \simeq 93.5\) MeV [90]. Similarly to the running coupling, accounting for higher orders in the running quark mass has no significant impact on our results.

3.1.2 Fit for \(3 \le X \le 6\)↩︎

The fitting functions provided in the following all depend on the quantity \(\mu_{B,0}(X)\) which represents the threshold for the vanishing of the strange quark density at different renormalization scales \(X\): \[\label{eq:muB0} \mu_{B,0}(X)=0.5565172 +\frac{0.7615706}{X^{1.0474315}} + 0.0044533\; X^{1.3455082} - 0.0074282\;X^{1.3726228}.\tag{16}\] It also gives the validity range: \(\mu_{B}\otimes X =\Lambda/\mu\in[\mu_{B,0},3.6] {\rm GeV} \otimes[3,6]\). For the threshold value \(\mu_{B}=\mu_{B,0}(X)\), the pressure and baryon density read \[\begin{align} P^{\rm th}(X)=&0.6607613- \frac{0.0016532}{X^{0.4460819}}-0.6595446 X^{0.0004152}, \\ \rho^{\rm th}_{B}(X)=&0.000092256-\frac{0.0041293}{X^{2.3736665}} + \frac{0.0013640}{X^{0.9529593}}. \end{align}\] To ensure the consistency of the thermodynamic relation between pressure and baryon density \(\partial_{\mu_{B}} P\equiv \rho_B\), both quantities have been fitted together, yielding \[\label{eq:Fit36} \begin{align} \frac{P^{\rm TC}_{\rm RGOPT}}{\rm GeV^4}\equiv\;&P_{\rm fit}(\tilde{\mu}_B=\frac{\mu_{B}}{\rm GeV},X)= P^{th}(X) + (\tilde{\mu}_B-\mu_{B,0}(X)) \rho_B^{th}(X)+(\tilde{\mu}_B -\mu_{B,0}(X))^2 \left(\frac{a_1\;X^{a_2}}{1 + a_3\;\tilde{\mu}_B^{a_4} X^{a_5} + a_6\; \tilde{\mu}_B^{a_7} X^{a_8}}\right)\\ &+ (\tilde{\mu}_B -\mu_{B,0}(X))^3 \left(\frac{a_9\;X^{a_{10}}}{1 + a_{11}\;\mu_{B}^{a_{12}} X^{a_{13}} + a_{14}\;\tilde{\mu}_B^{a_{15}} X^{a_{16}}}\right) + (\tilde{\mu}_B -\mu_{B,0}(X))^4 \left(\frac{a_{17}\;X^{a_{18}}}{1 + a_{19}\;\tilde{\mu}_B^{a_{20}} X^{a_{21}} + a_{22}\;\tilde{\mu}_B^{a_{23}} X^{a_{24}}}\right) ,\\ \end{align}\tag{17}\] with the coefficient values \(a_i\) given in Table [tbl:tabX36]. The baryonic density is easily obtained from Eq. (17 ) by taking the explicit derivative with respect to \(\tilde{\mu}_B\), all other parameters or functions being independent of \(\tilde{\mu}_B\).

The Root Mean Square Error (RMSE) for this fit is 0.36% for the pressure, while for the baryonic density, it is 0.32%.

Finally, the energy density is given by \[\frac{\mathcal{E}_{\rm RGOPT}}{\rm GeV^4}\equiv \mathcal{E}_{\rm fit}(\tilde{\mu}_B,X)=\;- P_{\rm fit}(\tilde{\mu}_B,X)+\mu_{B}\, \rho_B \,.\] We emphasize that in those fitted formulas all the highly nontrivial dependence on the running coupling \(\alpha_S(X \mu_B)\) and masses \(\tilde{m}_s(X \mu_B)\), \(\tilde{\eta}(X \mu_B)\), accounted for in our actual calculation through Eqs. (14 ), (15 ), (8 ), (7 ), is hidden and consistently embedded in the fitted \(X\)-dependence.

3.1.3 Fit for \(2\le X \le 3\)↩︎

To achieve sufficient accuracy, we considered an independent fit for lower \(2\le X\le 3\) values: similarly to the previous ones, the following fitting functions depend on \(\mu_{B,0}(X)\), the threshold for the vanishing of the strange quark density at different renormalization scales \(X\): \[\label{eq:muB0951} \mu_{B,0}(X)=0.0803569 + \frac{2.596390}{X^{5.125701}}-\frac{0.4668653}{X^{2.225740}} + \frac{1.1201412}{X^{0.380615}}\,.\tag{18}\]

The corresponding range of validity is now: \(\mu_{B}\otimes X =\Lambda/\mu\in[\mu_{B,0},3.6] {\rm GeV} \otimes[2,3]\). For the threshold value \(\mu_{B}=\mu_{B,0}(X)\), the pressure and baryon density read \[\begin{align} P^{\rm th}(X)=&0.0001877 - \frac{0.0205381}{X^{3.373904}} \; + \frac{0.8637938}{X^{2}} - {0.8505441}{X^{1.988330}} \\ \rho^{\rm th}_{B}(X)=& \theta(X-2.25)\left(0.2599325-\frac{0.0112160}{X^{3.653219}} -0.2590002\, X^{0.0026900} + 0.000108\, X\right)\,, \end{align}\]

where \(\theta(X)\) is the Heaviside function. Again, for thermodynamic consistency to hold between the pressure and baryon density, both quantities were fitted together, reading \[\begin{align} \frac{P^{\rm TC}_{\rm RGOPT}}{\rm GeV^4}\equiv\;&P_{\rm fit}(\tilde{\mu}_B=\frac{\mu_{B}}{\rm GeV},X)= P^{th}(X) + (\tilde{\mu}_B-\mu_{B,0}(X)) \rho_B^{th}(X)+(\tilde{\mu}_B -\mu_{B,0}(X))^{2.06306} \left(\frac{a_1\;X^{a_2}}{1 + a_3\;\tilde{\mu}_B^{a_4} X^{a_5}}\right)\\ &+ (\tilde{\mu}_B -\mu_{B,0}(X))^{2.940057} \left(\frac{a_6\;X^{a_{7}}}{1 + a_{8}\;\mu_{B}^{a_{9}} X^{a_{10}}}\right) + (\tilde{\mu}_B -\mu_{B,0}(X))^4 \left(\frac{a_{11}\;X^{a_{12}}}{1 + a_{13}\;\tilde{\mu}_B^{a_{14}} X^{a_{15}}}\right)\,.\\ \end{align}\label{Eq:pocket95X952953}\tag{19}\] with \(a_i\) coefficient values given in Table [tbl:tabX23]. The Root Mean Square Error (RMSE) for this fit is 1.01% for the pressure, while for the baryonic density it is 1.50%.

If needed, the above pocket formulas may be compared (or normalized to) the standard expressions for a gas of non-interacting massless quarks given in the standard Fermi-Dirac limit: \[P_{\rm FD} = \frac{N_c N_f}{12 \pi^2}\left ( \frac{\mu_B}{3} \right )^4 \,,\] and \[\rho_{\rm FD} =\frac{N_c N_f}{3 \pi^2}\left ( \frac{\mu_B}{3} \right )^3\,.\]

4 Numerical Results↩︎

In this Section, we present the numerical results for quark and hybrid stars. To describe quark matter, we compare two approaches, the conventional pQCD and the RGOPT resummation, both at NNLO, and to describe hadronic matter we use relativistic mean field (RMF) models with the DD2 parameterization of Ref. [81], which produces a stiff EoS, and two soft EoSs: the SFHo parameterization of Ref. [79] and the DIDY parametrization which includes hyperons and isospin-dependent couplings [80].

The RGOPT results are obtained using the fitting formulas of Section 3, while for pQCD we use the fitting function of Ref. [14], which reproduces the thermodynamically consistent pressure of Ref. [13].

4.1 Quark stars↩︎

Here, we present results for strange (\(N_f=2+1\)) QSs, comparing the RGOPT and pQCD predictions. In Fig. 3, we show the normalized pressure as a function of the baryon chemical potential for \(2\le X \le 4\). When comparing the two approximations for the same \(X\)-range, it is evident that the NNLO RGOPT exhibits substantially reduced scale dependence compared to NNLO pQCD. Notice that both methods yield very similar results at high \(\mu_B\) values.

Figure 4 shows the normalized pressure as a function of the baryon chemical potential and the pressure as a function of the energy density for the values of \(X\) that reproduce maximum QS masses, \(M_{\rm max}\,=2\), \(2.3\) and \(2.6\, M_\odot\). The dots correspond to the values at the center for the maximum QS mass configuration. The values of \(X\) for each approximation, together with several QS properties, are given in Table ¿tbl:tab:tab95quarks?.

The first point to note from the table is that the RGOPT values of the scale parameter \(X\) are higher than those obtained with pQCD. At the same time, the values of the couplings on the surface \(\alpha_s^{\rm surf}\), predicted by both approximations, are comparable and lie near the upper limit for which perturbative calculations are generally considered reliable, namely \(\alpha_s(\Lambda\sim1 \rm{GeV})\sim 0.423\) [92].

Regarding the scale dependence, it is interesting to note, from Table ¿tbl:tab:tab95quarks?, that for the \(\Delta M_{\rm max}=(2.60-2.00) M_\odot\) range considered here, the maximum masses have an almost linear dependence on \(X\) with angular coefficient \(\Delta M_{\rm max}/\Delta X \simeq 0.588 \,M_\odot\) per unit of \(X\) for RGOPT and \(\simeq 0.983 \,M_\odot\) per unit of \(X\) for pQCD, thus showing that the former is more stable to scale variations. It is also worth pointing out that the compactness \({\cal C}_{\rm max} = M_{\rm max}/R_{\rm max}\) is approximately constant in both cases, namely, \({\cal C}_{\rm max} \simeq 0.24\) for the RGOPT and \({\cal C}_{\rm max} \simeq 0.25\) for pQCD ⁵. This result could be anticipated by noticing that, to reproduce such maximum masses, the RGOPT requires larger values of \(X\) and coupling values that are slightly lower than those of pQCD, as the table shows (see, e.g., \(\alpha_s^{\rm surf}\)).

Finally, it is equally important to note that at NNLO both approximations are consistent with the Bodmer-Witten hypothesis for stable strange matter [93]–[95], since the energy per baryon of bulk strange quark matter at zero pressure is found to be below \(930\) MeV.

In Fig. 5 we show the mass-radius relations obtained by solving the Tolman-Oppenheimer-Volkoff equations [98], [99], using the EoSs of Fig. 4. The figure also displays observational constraints on the masses and radii of pulsars PSR J0740+6620 [68], [70], PSR J0614-3329 [73], PSR J0437-4715 [72] and the most recent analysis of NICER data for pulsars PSR J0030+0451 [96]. Also shown are the compact object HESS J1731-374 [74] and the band corresponding to the low mass component of the event GW190814 [97]. In general, both approaches respect the observational constraints, with pQCD predicting the most compact QSs. Considering the mass-radius curves that produce \(M_{\rm max}=2M_\odot\), the two methods predict results which lie within the overlap region of the \(\approx 1.4M_{\odot}\) pulsars PSR J0614-3329, PSR J0437-4715 and PSR J0030+0451. For the curves yielding \(M_{\rm max}=2.3M_\odot\), all results fall outside the constraint from pulsar PSR J0614-3329. However, the pQCD curve remains within the overlap region of the PSR J0614-3329, and PSR J0437-4715 constraints, while the RGOPT curve is consistent only with the \(95\%\) confidence region of PSR J0030+0451. Finally, for \(M_{\rm max}=2.6\, M_\odot\), the pQCD EoS still predicts configurations consistent with the constraints of HESS J1731-374 and PSR J0030+0451. In contrast, the RGOPT predictions lie outside both of these constraints.

4.2 Hybrid stars↩︎

In this Section, we discuss some properties of hybrid stars with a quark core described by the RGOPT, comparing our results with predictions from pQCD. To investigate the hadron-quark phase transition, a Maxwell construction is applied. For the hadron phase, three different EoSs have been chosen: a soft nucleonic EoS (SFHo [79]), a soft hyperonic EoS (DIDY [80]) and a stiff nucleonic EoS (DD2 [81]). These choices span approximately the range of phase space that is consistent with the presently known constraints on the nuclear EoS, and also consider the possible onset of hyperons.

For each hadronic model, we build two EoSs: i) the star for which quark matter has started to nucleate at its center, corresponding to the smallest scale \(X\) considered. These stellar objects define the most massive hybrid stars that already contain a quark core, although quite small; ii) the hybrid star that touches the 95% (full) confidence interval of the 2-D distribution of the PSR J0740+6620 [68], [70]. Such a condition defines the largest renormalization scale for each model, giving rise to stars with the largest quark core that are still consistent with observations. To define these two scenarios, the scale parameter \(X\) for RGOPT(pQCD) increases from 1.99(2.17) to 2.59(2.80) for DD2, from 2.28(2.25) to 2.98(2.86) for SFHo, and from 2.01(2.00) to 2.88(2.82) for DIDY (see also Table ¿tbl:tab:hybrid?). The two limiting scale values define the bands shown in Fig. 6, where the normalized pressure, \(P/P_{\rm FD}\), for RGOPT and pQCD is plotted as a function of the baryon chemical potential. The highest pressure corresponds to the largest scale parameter value.

Despite being an approximation, the Maxwell construction used to build the hadron-quark phase transition is known to produce reasonable results if the surface tension of hadronic matter in quark matter is large [100] although the actual value of the surface tension remains unsettled, see also [101]. We have plotted the hybrid EoS and the corresponding speed of sound squared in Fig. 7 for the three hadronic models, respectively, in the top and bottom panels. In this Figure, the dots identify the values at the center of the maximum mass star configuration. For DD2, the density at the center lies in the range \(\sim\) \(5.6-6.3\, \rho_0\), while for SFHo and DIDY the range spans \(6.4-8.3\, \rho_0\). The phase transition for the largest \(X\) occurs at \(\sim 2\, \rho_0\) for DD2, and \(\sim 4\, \rho_0\) for the other two models (thin lines). The strength of the phase transition, defined by the magnitude of the energy density jump, is quite weak. Interestingly, the hybrid EoSs constructed by matching the SFHo and DIDY hadronic models to either RGOPT or pQCD are nearly identical for the largest \(X\) values. For the smallest scale parameter value (thick lines), the phase transition occurs at densities of the order of \(4-6\, \rho_0\), the smallest value corresponding to DD2. This feature, which can be clearly seen in the bottom panels of Fig. 7, corresponds to the densities where the speed of sound jumps from the hadron branch to the quark branch. At these particular scales, the baryon density jump between the two phases is large, of the order or larger than \(2\, \rho_0\). Notice that, for DIDY, the phase transition with the largest \(X\) value occurs before the hyperon onset in the hadron branch, while for the smallest \(X\) value, the transition occurs after the hyperon onset. One should remark that it is the hyperon onset that defines the bump in the DIDY plot of the speed of sound squared (bottom right panel). Note also that, at the hadron-quark phase transition, the coupling constant \(\alpha_s^t\) is below or at the limit where perturbative expansions are considered valid (see Table ¿tbl:tab:hybrid?). For the quark phase, the pocket formula of Ref. [14] (pQCD) and the one derived in Sec. 3.1 for RGOPT were used to generate the corresponding plots.

The mass-radius curves for hybrid stars, together with the observational constraints from NICER and the compact object HESS J1731 − 347 are displayed in Fig. 8, while some of the main associated properties are summarized in Table ¿tbl:tab:hybrid?. The maximum mass configurations, corresponding to the smallest \(X\) value and the smallest quark core, are essentially defined by the properties of the hadronic EoS (namely SFHo, DIDY or DD2). Although the curves present a cusp like behavior, they contain a finite but small core of quark matter. In contrast, the configurations with the highest \(X\), correspond to the largest quark core predicted by each specific pair of EoS. Interestingly, the EoS with the largest quark core, corresponding to DD2+pQCD is in agreement with the observational constraint given by HESS J1731 - 347, as can be seen in the left panel of Fig. 8. As pointed out previously, for the hybrid-star configurations constructed using the SFHo and the DIDY EoS, the RGOPT and pQCD descriptions of quark matter produce nearly identical mass-radius relations.

In Fig. 9, the density profiles of the stars defined by the EoS with the highest \(X\) values are represented as functions of the distance to the surface. The star corresponding to the highest scale contains a core of quarks with a thickness of the order of 8 km (DD2) and 5 km (SFHo) or slightly smaller (DIDY), see Table ¿tbl:tab:hybrid?. The difference between DD2 and the other two models is easily understood since a stiff EoS favors a transition to quark matter at lower densities, allowing a large quark core to be more easily supported.

We have also calculated the density profiles of stars described by the quark pQCD EoS and the same models for hadron matter. The properties of these stars are similar to those obtained with the RGOPT EoS, with the main differences occurring if the hadron matter is described with the DD2 model. In this case, the \(X\) parameter is about 0.15 to 0.20 larger for the pQCD quark cores, and the star with the largest quark core compatible with the observation of the pulsar J0740 only has a small crust of hadron matter with a thickness of the order of 0.5 km. The onset of quark matter at a very low density has a direct impact on the radius of low mass stars, and the radius of a 1.4\(M_\odot\) star is about 1 km shorter for stars with a pQCD core.

5 Conclusions↩︎

We have extended a recent renormalization-group-induced resummation of pQCD, RGOPT at NNLO [65], to incorporate \(\beta\)-equilibrium and charge neutrality in a thermodynamically consistent EoS, tailored to describe massive quarks in the \(N_f =2+1\) case. Having such an EoS, from a first principles QCD evaluation, we have investigated the properties of (pure) QSs as well as the quark core of hybrid stars. Maximum masses corresponding to \(2, 2.3\) and 2.6 solar masses have been considered as phenomenological inputs to possibly constrain the corresponding arbitrary renormalization scale. In the case of pure QSs, our results show that the NNLO RGOPT predictions are less sensitive to scale variations than those generated by pQCD evaluations at the same perturbative order. The resummed RGOPT EoS final expression being involved, we have provided a compact formula that can be readily used in applications aiming to describe the quark sector associated with compact stellar objects simply by using \(\mu_B\) and \(X\) as inputs. As our results indicate, this formula incorporates all important RG properties so that the uncertainties (due to scale dependence) observed within its pQCD counterpart can be mitigated.

Our numerical analysis started by considering the case of pure QSs. The renormalization scale parameter was chosen by enforcing the corresponding QSs to be compatible with astrophysical observations and, as a consequence, values in the range \(X=3.08-3.58\) were obtained. It should be remarked that a larger value (\(X=4.01\)) needs to be considered if the low mass object of the observation GW190814 is a NS. At the same time, the range of possible \(X\) values obtained with a pQCD quark EoS, upon imposing the same observations, turned out to be about 40% smaller, 2.95 to 3.26 (and 3.56, if GW190814 is considered), as a direct consequence of a larger scale dependence. In addition, the lowest \(X\) value is \(\sim 5\%\) smaller.

To describe hybrid stars, we have used the RGOPT or pQCD EoS for the quark phase and merged it with three different EoSs that describe the hadron phase with different levels of stiffness (SFHo, DIDY, and DD2) and hadron content (since DIDY includes hyperons, in contrast to SFHo and DD2). To take into account the hadron-quark phase transition, the standard Maxwell construction has been adopted. Then, for each model, we generated two hybrid EoS extremes by varying \(X\): i) the lowest value describes a hybrid star with a small quark core and gives the most massive star, which is essentially defined by the hadron model; ii) the other extreme, which defines the largest \(X\), describes a hybrid star with the largest quark core, still compatible with the data of PSR J0740+6620. These two extremes set the range of interest for the scale parameter:

DD2 with RGOPT (pQCD), gives \(X=2.59-1.99\) (\(X=2.80-2.17\)) for maximum masses in the range \(1.84-2.29\, M_\odot\) (1.86 and \(2.23\, M_\odot\)). In the case of the hadronic SFHo and DIDY EoSs, together with the quark RGOPT or pQCD EoS, similar results were obtained: for SFHo, the maximum stellar masses lie within the range \(1.81-2.00\, M_\odot\) (\(X=2.98-2.28\)) and from \(1.82\) to \(2.15\, M_\odot\) (\(X=2.88-2.01\)) for DIDY.

Considering pQCD+DD2, hybrid stars with a small hadronic crust about 1 km thick were obtained. These stars are compatible with the compact object HESS J1731-347; similar conclusions were drawn in [102], [103]. Otherwise, stars with the largest quark core compatible with PSR J0740+6620 have a hadronic shell that is about 3 to 5 km thick and are not compatible with HESS J1731-347. The highest value of the speed of sound occurs in the hadronic branch and goes above the conformal limit, \(\sqrt{1/3}\). After the hadron-quark phase transition, the speed of sound drops to values below \(\sqrt{1/3}\) in the quark phase.

In the future, the possible existence of colorsuperconducting phases will be investigated. It will also be interesting to analyze the consequences of constraining the NS EoS with a smaller \(X\) range as the one obtained in the present study, when thermodynamical constraints, as described in [30], are imposed.

6 List of contributions to the \(N_f=2^*+1^*\) RGOPT pressure↩︎

For completeness, rather than repeating many long expressions, we simply indicate precisely where to find the relevant contributions of the complete NNLO pressure of Eq. (1 ) in [65], where the equation numbers refer to the later reference:

• \(P^{m}_{\rm LO, NLO}(m_i,\mu_i)\), \(P^{v}_{\rm LO, NLO}(m_i)\), respectively in Eqs. (2), (3), (23), (24);

• \(P_{\rm 2GI}^{N_f=2^*+1^*}(m_i,\mu_i)\), \(P_{\rm VM}^{N_f=2^*+1^*}(m_i,\mu_i)\) in Eqs. (8), (9), (18), (46),(47);

• \(P_{\rm Ring}^{N_f=2^*+1^*}(m)\equiv -\Omega_{\rm Ring}^{N_f=2^*+1^*}(m)\) in Eq. (21);

• \(P^v_{\rm NNLO,i}(m_i)\) in Eqs. (25), (26), (44);

• Finally \(P_{sub,i }\), \(P_{sub}^{nd}\) in Eqs. (45), (A.9), (A.10).

7 Renormalization group material↩︎

Since it is a central quantity in our approach, we recall the expression of the massive (homogeneous) renormalization group operator \[\begin{align} \Lambda \frac{d}{d \Lambda}\equiv \Lambda \frac{\partial}{\partial \Lambda} + \beta(g_s^2 )\frac{\partial}{\partial g_s^2 } - \gamma_m(g_s^2 ) m \frac{\partial}{\partial m} \; \label{RG} \end{align}\tag{20}\] where in our NNLO analysis we use the beta-function \(\beta(g_s^2)\) and quark mass anomalous dimension gamma-function \(\gamma_m(g_s^2)\) at three-loop order, defined in our conventions as \[\beta\left(g_s^2 \right)=-2g_s^4\left(b_0 +b_1g_s^2+b_2 g_s^4+ \mathcal{O}\left(g^6\right) \right)\;,\] \[\gamma_m\left(g_s^2\right)=g_s^2\left( \gamma_0+\gamma_1g_s^2 +\gamma_2 g_s^4 + \mathcal{O}\left(g^6\right)\right)\;,\] with the coefficients, for the relevant QCD case with \(C_A=N_c\), \(C_F=4/3\), \(N_c=3\), and \(N_f\) quark flavors: \[\label{eq:BetaCoef} \begin{align} b_0= & \frac{1}{(4\pi)^2}\left({\frac{11}{3}\,}C_A-\frac{2}{3}N_f\right) \;\;,\;\; \end{align}\tag{21}\] \[\label{eq:GammaCoef} \begin{align} \gamma_0 = & \frac{2}{(4\pi)^2} (N_c C_F ), \;\;\; \end{align}\tag{22}\] and higher order QCD RG coefficients up to relevant three-loop order given in our normalization in [65] (see e.g. refs[104]–[106] for QCD RG coefficients up to five-loop order). In addition, as mentioned in Sec. II, upon including massive vacuum terms the renormalization group invariance of the massive pressure requires additional zero-point subtraction contributions, \(P_{sub}\sim \sum_k s_k (g_s^2)^{(k-1)}\) in Eq.(1 ), directly related to the fact that RG invariance actually involves [52], [82], [83] the vacuum energy anomalous dimension \(\hat{\Gamma}^0(g_s^2)\), i.e.: \[\Lambda \frac{d}{d \Lambda} P(m,g_s^2) -m^4 \hat{\Gamma}^0(g_s^2) \equiv 0\] with \[\hat{\Gamma}^0(g_s^2) \equiv \sum_k \Gamma^0_k (g_s^2)^{k}\] In the normalization of the pressure in Eq. (1 ), the resulting subtraction coefficients \(s_k\) are \[\label{eq:s95i32Quark} \begin{align} s_0 \equiv &\frac{\Gamma^0_0}{2\big(b_0-2 \gamma_0\big)}=\frac{3}{7} \;\;, \;\; \end{align}\tag{23}\] with higher order coefficients given in [65].
Finally, one should bear in mind that performing a renormalization scheme change (RSC) according to Eq. (6 ) also implies consistent \(B_2\)-dependent modifications [65] in the higher order RG coefficients \(\gamma_k, s_k\) above.

References↩︎