Recovering the conformal limit of color superconducting quark matter within a confining density functional approach

We generalize a recently proposed confining relativistic density-functional approach to the case of density dependent vector and diquark couplings. The particular behavior of these couplings is motivated by the non-perturbative gluon exchange in dense quark matter and provides the conformal limit at asymptotically high densities. We demonstrate that this feature of the quark matter EoS is consistent with a significant stiffness in the density range typical for the interiors of neutron stars. In order to model these astrophysical objects we construct a family of hybrid quark-hadron EoSs of cold stellar matter. We also confront our approach with the observational constraints on the mass-radius relation of neutron stars and their tidal deformabilities and argue in favor of a quark matter onset at masses below $1.0~{\rm M}_\odot$


I. INTRODUCTION
Modern multi-messenger observations of neutron stars (NSs) and their mergers provide new measurements of their masses and radii. These data are important constraints on the equation of state (EoS) of cold, dense matter in a region of the QCD phase diagram which is inaccessible to ab initio simulations of QCD on the lattice or heavy-ion collision experiments. The results of the analysis give at a mass of 1.4 M a radius R 1.4 = 11.7 +0.86 −0.81 km [1] and at 2.0 M a radius of R 2.0 = 13.7 +2.6 −1.5 km [2]. These results imply that the neutron star matter EoS should not be too stiff at densities below twice the saturation density n 0 = 0.15 fm −3 (roughly corresponding to the central density of a NS with M = 1.4 M ), but has to be stiff enough at higher densities to allow for a maximum mass above 2.0 M . These new constraints on the NS mass-radius relation can be fulfilled within the purely nucleonic scenario for the NS interiors [3]. At the same time, approaches based on realistic nuclear interactions imply appearance of hyperons in the NS interiors, which softens EoS of nuclear matter and lowers the NS maximum mass M max . For example, ab initio Brueckner-Hartree-Fock and cluster variational methods with the microscopic interaction potentials fitted to the nucleusnucleus scattering data and properties of hypernuclei yield M max barely reaching 2.0 M [4]. The analysis performed within a set of EoSs derived from relativistic density functional theory constrained by the results of chiral effective field theory, terrestrial experiments, astrophysical observations and reproducing hyperon potentials in the symmetric nuclear matter at saturation density provides marginal agreement with the NICER constraints on the NS maximum mass [2,5]. Stiffer hadronic EoSs provide better agreement even in the presence of hyperons, e.g. DD2 EoS with hyperons yields M max = 2.1 M [6]. However, stiff hadronic EoSs are discriminated by the requirement that tidal deformability of a 1.4 M NS falls within the range Λ 1.4 = 70−580 extracted from the analysis of GW170817 [7]. These complications are naturally removed when the scenario with a low onset density for the transition to stiff quark matter in the NS core is considered, so that all the above conditions can be fulfilled simultaneously. It is important to note, that recent model agnostic statistical analyses report the viability of EoSs without a strong first order phase transition [8,9]. In this case the sharp interface between quark and hadron matter is smoothened, e.g. by inhomogeneous pasta structures [10]. Moreover, an early onset of deconfinement for star masses around 0.5 M could at present not discriminated observationally from a sequence without a phase transition, even with the recent measurement of a strangely light neutron star [11]. This spectacular measurement is in excellent agreement with the scenario of an early onset of deconfinement.
As it has been discussed in detail in the recent review by Baym et al. [12], a NS EoS with stiff quark matter requires a repulsive vector meanfield and strong color superconductivity with sufficiently large diquark pairing gap for the early onset of the deconfinement transition. When one aims at a sufficiently general formulation of the quark matter EoS which should also be suitable for the description of systems at finite temperatures like in supernova explosions or neutron star mergers, a confining relativistic density functional approach has proven successful [13]. The model developed in [13] has recently been generalized in [14] so that its Lagrangian obeys chiral symmetry and describes color superconductivity.
However, in these approaches, the vector meanfield persists at high densities and thus the quark matter EoS remains stiff with a squared sound speed well exceeding the conformal limit value of c 2 s = 1/3. Many authors do not recognize this situation as a problem since the densities at which perturbative QCD (pQCD) provides a reliable EoS model are about one order of magnitude larger than the central densities in the most massive NSs. Nevertheless, it has been shown recently [15] that pQCD actually can constrain the EoS at NS densities, just by demanding thermodynamic stability and causality. Therefore, in the present work we want to present a possible generalization of the RDF approach to dense quark matter which recovers the conformal limit at high densities.
The organization of the paper is as follows. In the next section we generalize the RDF approach developed in Ref. [14] to the case of density dependent vector and diquark couplings. The EoS of cold, color-superconducting quark matter is modelled in Section III. Its convergence to the conformal limit is analyzed in the same section. Section IV is devoted to application of the developed EoS to modelling compact stars with quark cores. The results are summarized and discussed in Section V.

II. CONFINING RDF APPROACH WITH DENSITY DEPENDENT VECTOR AND DIQUARK COUPLINGS
A generalization of the RDF approach for the description of color-superconducting quark matter to the case of density dependent vector and diquark couplings can be performed within the Lagrangian formalism developed in Refs. [14,16]. In the case of two quark flavors the fundamental dynamical variables of the approach are quark fields represented by the flavor spinor q T = (u, d). We note that a first application of the present approach to the three-flavor case was given in Ref. [17]. The interaction terms are chosen in the contact current-current form (qΓq) 2 withΓ being an interaction vertex. In the case of scalar (Γ = 1) and pseudoscalar (Γ = iγ 5 τ) channels the corresponding quark bilinears are composed to the chirally symmetric combination (qq) 2 + (qiγ 5 τq) 2 providing the corresponding symmetry of interaction. The model Lagrangian can be written as where m is the current mass of the two light quark flavors. The potential U accounts for the attractive chirally symmetric interaction in scalar and pseudoscalar channels where constants D 0 and α control the interaction strength and constituent quark mass in the vacuum [14,16], respectively, while qq 0 is the vacuum value of the chiral condensate. In what follows the subscript index "0" denotes the quantities defined in the vacuum. The present parameterization of U is motivated by the string-flip model (SFM) [18,19], which assumes that the interparticle interaction energy is proportional to mean separation between quarks. The model Lagrangian includes terms representing vector repulsion and diquark pairing interactions where charge conjugated quark field is q c = iγ 2 γ 0 q T and A = 2, 5, 7 labels the antisymmetric generators λ A /2 of the SU(3) color group, so that the ansatz (4) for the diquark current fulfills the Pauli principle for the quark pair. These interaction channels are important for compact star phenomenology [12]. In Refs. [14,16,20] the couplings G V and G D were set to be constants. Here we consider them as medium dependent functions. This section presents a general treatment, while the specific parameterization of the vector and diquark couplings adopted in this work is considered in Section III. The naive introduction of a medium dependence for G V and G D can break thermodynamic consistency by violating thermodynamic identities similarly to the case of the naive introduction of a medium dependent dispersion relation [21]. In order to circumvent this problem we follow the strategy of Ref. [22] and introduce the so called rearrangement terms Θ V and Θ D to Eqs. (3) and (4). Similar to G V and G D , they are some medium dependent functions which should be defined in agreement with the corresponding couplings. These rearrangement terms vanish at constant vector and diquark couplings. Their signs in Eqs. (3) and (4) are conventional. The present choice is motivated by the fact that the corresponding terms in the Lagrangian represent repulsive and attractive interactions.
Expanding the potential U around the mean-field expectation values qq 0 and qiγ 5 τq = 0 up to the second order terms and inserting the result to Eq. (1) yields an effective Lagrangian. At this order the only non-vanishing expansion coefficients are Hereafter the subscript index "MF" denotes the quantities defined by the mean field approximation. The resulting effective Lagrangian has the current-current interaction form of the NJL type models, where m * = m + Σ MF is the constituent quark mass. It follows from this effective Lagrangian that Σ MF is nothing else than scalar self-energy of the quarks at the mean-field level, while G S and G PS correspond to the effective couplings of quark interaction in the scalar and pseudoscalar channels, respectively. These couplings do not coincide in the general case. This corresponds to an explicit violation of chiral symmetry which results from expanding the Lagrangian L around the mean-field solution, which is known to be chirally broken. At the same time, the dynamical restoration of chiral symmetry at high temperatures and densities leads to the asymptotic coincidence of G S and G PS [14,16]. With the effective Lagrangian L eff , the partition function can be represented as a functional integral over quark fields where integration over the Euclidean space-time is limited to the inverse temperature 1/T ≡ β =´dτ and the volume V =´dx.
The diagonal matrixμ = diag(µ u , µ d ) stands for the quark chemical potentials. They can be expressed through the baryonic µ B and electric µ Q chemical potentials as where subscript index f = u, d labels quark flavors and Q f is their electric charge. The next step corresponds to bosonizing the partition function by means of the Hubbard-Stratonovich transformation. This introduces collective scalar (σ), pseudoscalar ( π), vector (ω µ ) and complex scalar diquark (∆ A ) fields. They are coupled to the corresponding bilinears of quark fields, qq − qq , qiγ 5 τq, qγ µ q and qiγ 5 τ 2 λ A q, respectively. It is worth noticing that the medium dependence of the couplings G S , G PS , G V and G D does not affect this procedure since none of them includes any dynamical variable. It is convenient to treat the bosonized partition function within the Nambu-Gorkov formalism. Here we just outline the main aspects of the formalism and summarize the results. The interested readers are referred to Refs. [14,23]. In this case quark fields are collected to the Nambu-Gorkov bispinor Q T = (q q c )/ √ 2, while the partition function becomes Here the propagator of the Nambu-Gorkov bispinors reads with S −1 ± = i/ ∂ ± / ω − m * ± γ 0μ . The exponential in the second line of Eq. (10) is nothing else than the quark-meson part of the bosonized action. The quark fields enter this action quadratically and can therefore be integrated out analytically yielding Tr ln βS −1 /2 in the exponential. The trace Tr hereafter is performed over the color, flavor, Dirac, three-momentum and Matsubara indices. The last ones appear after going over to the momentum representation which yields S −1 ± = / k − m * with k 0 = iz n ±μ * , z n = (2n + 1)πT defining a fermionic Matsubara frequency and µ * f = µ f + ω being effective chemical potential of quarks.
The action in Eq. (10) gives direct access to the Euler-Lagrange equations of the scalar, pseudoscalar, vector and diquark fields. Averaging these equations for the vector and diquark fields one obtains ω µ = −2G V qγ µ q and ∆ A = 2G D q c iτ 2 γ 5 λ A q , respectively. By a proper Lorentz transform the vector field average attains the form ω µ = g µ0 ω with ω = −2G V q + q . Furthermore, there exists a global color rotation which leaves ∆ 2 as the only diquark field with a nonvanishing expectation value at the mean field level. We note that only its modulus ∆ = |∆ 2 | appears in the expression for thermodynamic potential. Averaging the Euler-Lagrange equations for scalar and pseudoscalar fields shows that σ and π vanish at mean field [14]. Thus, σ and π have beyond the mean-field nature and represent the corresponding mesonic correlations of quarks. Within the Gaussian approximation the back-reaction of these correlations on the quark propagator is neglected [23]. This allows to expand Tr ln βS −1 up to the second order in σ and π. The second order terms are quadratic in the mean-field quark propagator S MF and thus represent one-loop polarization operators of (pseudo)scalar mesons. The latter can be used in order to construct mesonic propagators and to extract the corresponding masses from the position of the propagator poles. Within the generalized Beth-Uhlenbeck approach, mesonic propagators also can be used in order to obtain beyond mean-field contributions to the thermodynamic potential [16,23]. In the present work, however, they are neglected because we restrict ourselves to the meanfield approximation. For this we replace scalar, pseudoscalar, vector and diquark fields in Eq. (10) by their expectation values and reduce the corresponding functional integrals. This yields The first term in this expression is due to the contribution of quark quasiparticles It includes the single particle energies shifted by the effective chemical potential and distribution functions, i.e.
Here k f = √ k 2 + m * 2 and the subscript index c = r, g, b labels quark color states. The color vector ∆ c = (∆, ∆, 0) is introduced in order to unify the notations and a = ± distinguishes particles and antiparticles. The dispersion relation (14) can be obtained by solving det(S −1 MF ) = 0 with respect to the zeroth component of quark four momentum k. It shows that only red and green quarks are paired exhibiting the gap ∆ in their one-particle energy spectrum, while blue quarks are unpaired. The zero point terms in the expression for Ω q are regularized by smooth cut-off in the Gaussian form In Ref. [14] such a form was chosen in order to prevent a discontinuous behavior of various thermodynamic quantities which would have been obtained for the 2SC phase of quark matter with a sharp cutoff as g k = θ(Λ − |k|).
The thermodynamic definition of the number density of a given quark flavor corresponds to the thermodynamic identity f + f = −∂Ω/∂µ f , which should be used carefully since vector and diquark couplings are medium dependent functions. On the other hand, the statistical definition of this quantity implies f + f = −∂Ω q /∂µ f . The thermodynamic consistency of the present approach is provided when these two definitions coincide. Thus, we require where the mean-field equations ω = −2G V q + q and ∆ = 2G D | q c iτ 2 γ 5 λ 2 q | were used on the second step. Fulfilment of this condition requires the rearrangement terms Θ V and Θ D to be defined in accordance with the couplings G V and G D . The corresponding relations can be easily found by assuming that Θ V , G V , Θ D and G D are functions of q + q and | q c iτ 2 γ 5 λ 2 q |, respectively. In this case Eq. (16) leads to From these relations it is seen that the rearrangement terms, indeed, vanish if the couplings are constant. Using these relations number density of a given quark flavor, chiral condensate and modulus of the diquark one can be found from the quark part of the thermodynamic potential as We note that the Dirac delta-function in Eqs. (18) and (19) appears due to differentiating the sign-function from the dispersion relation (14). It is also worth mentioning that the definitions of the rearrangement terms given by Eq. (17) along with the mean-field equations ω = −2G V q + q and ∆ = 2G D | q c iτ 2 γ 5 λ 2 q | are sufficient in order to obtain the number density of a given quark flavor, the chiral condensate and the modulus of the diquark one in the form (18) - (20). This holds for any functional dependence of the couplings G V and G D on their arguments. With Eqs. (18) and (20) the mean-field equations for the vector field and diquark pairing gap can be given an explicit form. Furthermore, Eq. (19) should be understood as another mean-field equation with respect to chiral condensate. Its solution along with the solutions of the mean-field equations for vector field and pairing gap minimize the thermodynamic potential. For the readers convenience in Appendix A we explicitly analyze the conditions providing the minimum of Ω and derive from them the mean-field equations mentioned above as well as Eqs. (18) - (20). Once these mean-field equations are consistently solved, pressure, entropy and energy density can be found using the thermodynamic identities p = Ω 0 −Ω, s = ∂p/∂T and ε = f µ f f + f +T s− p, while squared speed of sound is defined as the derivative c 2 S = d p/dε calculated at constant entropy. Below we also analyse the dimensionless interaction measure δ = 1/3 − p/ε being nothing else than the trace of the energy momentum tensor scaled by the conformal limit for the pressure which is 3ε.
It is worth noticing that the present approach with density dependent vector and diquark couplings is equivalent to a density functional approach in the spirit of Ref. [13]. In the case of vector repulsion and diquark pairing the corresponding density functionals depend on the quark bilinears q + q and q c iτ 2 γ 5 λ 2 q, qiτ 2 γ 5 λ 2 q c , respectively. Consistency with the present approach with medium dependent vector and diquark couplings is provided by Expanding these potentials around the mean-field solutions up to the first order terms produces the vector and diquark selfenergies of quarks whereΣ D is an antidiagonal matrix in the Nambu-Gorkov space due to the fact that U D depends on two dynamical variables q c iτ 2 γ 5 λ 2 q and qiτ 2 γ 5 λ 2 q c . The vector self-energy shifts the quark chemical potential by exactly the same amount as ω = −2G V q + q . The corresponding pressure term −U V,MF + q + q Σ V also coincides with ω 2 /4G V + Θ V . The diquark self-energy coincides with the non-diagonal terms in the inverse Nambu-Gorkov propagator if the diquark fields in Eq. (11) are replaced by their expectation values discussed above. In this case the pressure term coming from the expansion of the diquark potential The present model has four parameters relevant to the QCD phenomenology, which are m, D 0 , α and Λ [14]. The pion mass M π and decay constant F π are the most important observables in this context. An analysis of the scalar mode mass M σ also was performed despite the fact that its experimental status is far from being clear. Our approach allows M σ in a wide interval covering the masses of all the experimental candidates. Typically, the lightest state f 0 (500) is considered as a candidate for the scalar meson role. It, however, has a large decay width of about 500-1000 MeV [24] and should be considered rather a tetraquark state than a traditional quark-antiquark meson [25]. Therefore, it is not appropriate to fit the vacuum parameters of the low-energy QCD model using the f 0 (500) state as a quark-antiquark meson. Our analysis uses instead the f 0 (980) state as a scalar meson. Our approach does not fit the vacuum value of chiral condensate per flavor | ll 1 GeV 0 | 1/3 = 241 MeV found from QCD sum rules at the renormalization scale 1 GeV [26]. This problem is typical for most of the chiral quark matter models [27]. Therefore, within the present approach we allowed the chiral condensate to have a somewhat larger value in order to have a reasonable value of the pseudocritical temperature T PC = 163 MeV defined by the peak position of chiral susceptibility. The model parameters defined using the above strategy along with the resulting physical quantities are presented in Table I. This parameter set yields m * = 718 MeV in the vacuum, which provides an efficient phenomenological confinement of quarks due to their high masses at low temperatures and densities.  The parameterization of the vector coupling adopted in this work is motivated by the analysis of the quark repulsion energy due to non-perturbative gluon exchange of QCD in the Landau gauge [28]. In the normal phase of symmetric quark matter it implies G V ∝ (9M 2 g + 8k 2 F ) −1 , with M g and k F = (π 2 q + q /2) 1/3 being the non-perturbative gluon mass and the quark Fermi momentum, respectively. With this we introduce At M g → ∞ this parameterization corresponds to constant vector and diquark couplings. At the same time, the solution of the gluon Schwinger-Dyson equations in the Landau gauge implies M g = 300 − 700 MeV [29,30]. The value of M g can be also estimated based on the Shifman, Vainshtein, and Zakharov expansion of the two-point current correlation functions within massive gauge invariant QCD [31]. For the frozen QCD structure constant α s = 0.2 and transferred momentum Q 2 = 10 GeV 2 that approach yields M g = 750 MeV. At α s = π, which is expected for the non-perturbative regime [32], the effective gluon mass becomes 516 MeV. At the same time, for the transferred momentum coinciding with the ultraviolet cut-off Λ from Table  I, i.e. for Q 2 = Λ 2 , one gets M g = 942 MeV at α s = 0.2 and M g = 792 MeV at α s = π. Below, we consider several values of the non-perturbative gluon mass, covering a range that contains the values mentioned above. This allows us to demonstrate continuous convergence to the conformal limit at all finite M g . We treat the vacuum values of the vector G V0 and diquark G D0 couplings as free parameters. They are parameterized as dimensionless ratios defined with the vacuum value of the scalar coupling G S 0 = 18.1 GeV −2 through η V ≡ G V0 /G S 0 and η D ≡ G D0 /G S 0 . In what follows, pairs of numbers (η V , η D ) are used in order to label different EoS parmetrizations obtained within the present model. It is necessary to stress that the physical values of η D are limited from above by the value η max , beyond which already the vacuum state would become color superconducting [14], see also [33,34]. For the chosen values of the model parameters η max D = 0.78. In order to not go to the marginal values of the diquark coupling we constrain the consideration to the range η D < 0.77.

III. COLD QUARK MATTER
Studying cold quark matter at vanishing temperatures is of practical interest for modeling compact stars. At T = 0 the single particle distribution functions of quarks reduce to unit step-functions, i.e. f + k f c = θ(− + k f c ) and f − k f c = 0. Therefore, the antiquark terms are absent at T = 0, while for the quark ones the integration over k is limited by the Fermi momentum defined by the condition + k f b = 0. Constructing the EoS of quark matter requires solving mean-field equations for chiral condensate, vector field and diquark pairing. Solutions of these equations give direct access to the single quark energies + k f c , which are needed in order to calculate the thermodynamic quantities mentioned in the previous section. Before considering these quantities in detail we would like to analyze the high density asymptotics of our approach and show that it is consistent with the conformal limit of strongly interacting matter. At high densities the zero point terms, the quark masses as well as all terms related to U can be neglected. In order to analyze the behavior of the pairing gap we notice that in the considered regime G D ∝ | q c iτ 2 γ 5 λ 2 q | −2/3 . Therefore, using the pairing gap equation ∆ = 2G D | q c iτ 2 γ 5 λ 2 q | and Eq. (20) we obtain where the summation over the color index was performed explicitly in the second step. For ∆ µ * f , the single particle energy of paired quarks can be approximated as + k f r −∆ and the bracket on the right hand side of this relation behaves as ∼ f µ * 3 f . This leads to ∆ ∼ ( f µ * 3 f ) 1/3 , which contradicts the original assumption. At ∆ ∼ µ * f or ∆ µ * f the bracket in Eq. (25) ∼ ∆ f µ * 2 f . Thus, the high density asymptotic of the pairing gap is ∆ ∼ ( f µ * 2 f ) 1/2 . The corresponding contribution to the pressure behaves as p D ≡ −∆ 2 /4G D − Θ D ∼ ∆ 4 . In order to show that p D scales as the vector field term p V ≡ ω 2 /4G V + Θ V ∼ q + q 4/3 we consider the quark number density The first term in the bracket of this expression contributes ∼ f ∆µ * 2 f to the quark number density. Following the above analysis of the pairing gap, this contribution is of the same order as the second term in the brackets in Eq. this we conclude that p D ∼ q + q 4/3 . Similarly, for the zero temperature quark pressure at high densities we obtain Thus, the total pressure at high densities is proportional to n 4/3 B with n B = q + q /3 being the baryon charge density. Along with the thermodynamic identity n B = ∂p/∂µ B this leads to the high density scaling p ∼ µ 4 B , which respects the conformal limit. Convergence to this limit is shown in Fig. 1. At small values of the baryonic chemical potentials chiral symmetry is broken, quark masses are high and the paring gap vanishes. The pressure also vanishes in this regime leading to a vanishing speed of sound and an interaction measure equal to one third. At a certain µ B , the quark mass discontinuously drops to some small value, while the pairing gap becomes finite. This discontinuous change signals a first order transition from chirally broken quark matter to the 2SC phase. At this transition the speed of sound jumps to some finite value, while the interaction measure exhibits a kink. In the general case, δ experiences a jump of the amplitude p c (1/ε 1 − 1/ε 2 ), where p c is pressure at the phase transition and ε 1 and ε 2 are energy densities of the coexisting low and high density phases. At zero temperature this jump of the interaction measure degenerates to kink with δ = 1/3 since p c = 0 in this case. Above the transition c 2 S decreases and δ has a minimum with some negative value. Position of this minimum can be found by requiring dδ/dε = p/ε 2 − c 2 S /ε = 0 or, equivalently, c 2 S = p/ε. This condition can be fulfilled only if p > 0 just because c 2 S are ε are positively defined. This explains why the minimum of δ is located in the region with positive pressure, i.e. above the transition from chirally broken quark matter to the 2SC phase. Qualitatively identical behavior of c 2 S and δ in the density range with positive pressure is exhibited by the parametric model of the quark matter EoS from Ref. [35]. At any finite value of the non-perturbative gluon mass and for µ B → ∞, the squared speed of sound and the interaction measure converge to their values in the conformal limit, c 2 S = 1/3 and δ = 0, respectively. It is important to stress, that vanishing of the interaction measure at finite µ B does not signal about reaching the conformal limit, which also requires ∂δ/∂µ B = 0. Note, these two conditions necessarily lead to c 2 S = 1/3 not being the case of finite values of the baryonic chemical potential. The smaller M g the faster is this convergence. In other words, the non-perturbative gluon mass defines the scale at which the quark-quark interaction effects cease. It is also worth mentioning that within the considered model c 2 S → 1/3 from above and δ → 0 from below. This means that vector repulsion between quarks dominates the attractive pairing interaction even at high densities. An alternative scenario with c 2 S approaching the conformal limit from below and δ approaching it from above requires domination of the attractive (pairing) interactions at high µ B . Such a case can be provided, e.g., if the common gluon mass parameter M g in Eqs. (23) and (24) is replaced by two independent vector (M gV ) and diquark (M gD ) masses with M gV < M gD . In this case the repulsive vector interaction ceases out before the attractive pairing one. An analysis of this scenario is beyond the scope of the present work. It is also important to note that for finite M g the variation of c 2 S in the density region typical for NSs is much more pronounced than for M g → ∞ which corresponds to the case of constant vector and diquark couplings considered in Ref. [14]. As is seen from Fig. 1, in symmetric quark matter this variation within the interval of µ B from 1 GeV to 2 GeV is about 10%.
Qualitatively the same behavior of c 2 S and δ is observed at any values of η V and η D and in the case of electrically neutral β-equilibrated quark matter, which is important for modelling NSs. Electric neutrality requires a proper amount of electrons with chemical potential µ e = µ u − µ d providing β-equilibrium. At small densities where quarks are confined we construct a phase transition of the quark matter EoS with the hadron one. We use DD2 EoS with hyperons, which is referred to as DD2npY-T [6]. The quark and hadron EoSs are matched by means of the Maxwell construction corresponding to the first order phase transition. Such a hybrid EoS is shown in Fig. 2. It is worth mentioning that the hadron-to-quark matter transition happens above the transition from chirally broken phase to the 2SC one in pure quark matter. Therefore, the quark matter branch of hybrid EoS is already color superconducting. For the considered values of M g diquark coupling strongly influences the onset density of quark matter decreasing with η D . The stiffness of the quark matter EoS is regulated by the vector coupling. A correlated variation of η V and η D allows us to generate a family of quark-hadron EoSs consistent with the constraints obtained within the multipolytrope analysis of the observational data of PSR J1614+2230 [36] and PSR J0740+6620 [2] as well as the statistical analysis from Ref. [37].
We want to point out that different values of the non-perturbative gluon mass lead qualitative different properties of the resulting EoS. At M g = 500 MeV the onset density of quark matter is limited from below by about 450 MeV fm −3 corresponding to η V = 0 and η D = 0.77. Smaller onset densities can be obtained only if going to the unphysical region η D > η max D . At physical values of the diquark coupling, the onset of quark matter for M g = 500 MeV always occurs after the hyperonization of the matter. As a result, the lower limit for the NS maximum mass can be reached only marginally, while the tidal deformability is outside the observational bounds. Already at M g = 600 MeV the onset density of quark matter is not limited from below, which provides a positive feedback to the problem of fulfilling the observational constraints. For the considered values of the vector and diquark couplings and the lower value of the gluon mass, M g = 500 MeV, the quark-hadron mixed phase is located at ε = 500 − 900 MeV fm −3 . For M g ≥ 600 MeV, however, the mixed phase lies at lower energy densities ε = 180 − 500 MeV fm −3 . As it is shown below, small values of the non-perturbative gluon mass are inconsistent with the observational constraints on the NS mass-radius relation and the bound on the tidal deformability. Based on this we conclude that our approach predicts the quark-hadron transition at energy densities below 500 MeV fm −3 . This range coincides with the lattice QCD results related to the chiral crossover region at vanishing chemical potential [38]. More recent analyses of the modern NS mass and radius constraints from multi-messenger observations using hybrid EoS with color superconducting quark matter do also find the hadron-to-quark matter transition at energy densities below 500 MeV fm −3 , see [39,40]. Such a coincidence of energy density domain for the phase transformation in two very different regions of the QCD phase diagram has already been observed in [41]. However, mechanism of this transformation at finite temperatures and finite densities can substantially differ. This makes a theoretical interpretation of the above "universality" of the transition energy density challenging.
It is worth mentioning that for M g ≥ 600 MeV our approach is able to generate EoSs, which at high densities are significantly stiffer than the hadronic EoS DD2npY-T. At the same time, this feature is consistent with approaching the conformal limit c 2 S → 1/3. Indeed, as is seen from Fig. 3, after the hadron-to-quark matter transition the squared speed of sound has a value of ∼ 0.5 and then decreases, approaching the conformal value at asymptotically high densities. Such a decrease of c 2 S in the density range typical for NSs was recently reported based on the model agnostic statistical analysis [42]. Within our approach the speed of sound reaches its maximal value at the quark boundary of the quark-hadron mixed phase. In other words, c 2 S peaks in the color superconducting 2SC phase right after the hadron-to-quark matter transition. For the EoSs providing the best agreement with the observational constraints discussed below (red and purple curves on the upper right and lower panels in Fig. 3), this maximum is located at ε 300 − 400 MeV fm −3 . Remarkably, this range of energy densities corresponding to the speed of sound maximum is in a very good agreement with the results of Ref. [43] obtained within the model agnostic statistical analysis, which also respects the conformal limit. A similar peak of c 2 S was also reported in Ref. [8] for the scenario of a first order phase transition FOPT-1 with the Group 1 of EoS.
For M g ≥ 600 MeV and energy densities interesting for the phenomenology of NSs, c 2 S varies around 0.45-0.55. Such values have been obtained for color superconducting quark matter within the nonlocal Nambu-Jona-Lasinio model with covariant [44] and also with instantaneous formfactors [40]. At the same time, for a given EoS of quark matter, the relative variation of c 2 S is about 10 %, being in tension with the assumption of the constant speed of sound (CSS) parameterization of the quark matter EoS [45,46].

IV. COMPACT STARS WITH QUARK CORES
Observational data on the masses and radii of NSs give important constraints on their EoS. These constraints include the measurement of the lower limit of the TOV maximum mass 2.01 +0.04 +0.04 M for the pulsar PSR J0348+0432 in a binary system with a white dwarf companion [47], the results of the Bayesian analysis of the observational data from PSR J0740+6620 [2,5] and PSR J0030+0451 [48,49], as well as the constraints on the masses and radii obtained from the gravitational wave signal and the kilonova light curve of the binary neutron star merger GW170817 [7,50,51]. We confronted these constraints with the massradius relations obtained by solving the problem of relativistic hydrostatic equilibrium, i.e. the TOV equation supplemented with the necessary boundary condition. The family of hybrid quark-hadron EoSs presented in Fig. 2 was used as an input for this task. The corresponding mass-radius relations are shown in Fig. 4. For M g = 500 MeV, the quark matter onset mass is rather high, while the TOV maximum masses reach the observational limit, if at all, only marginally. A similar problem arises when one compares the value for R 1.4 , the radius for a NS with a mass of 1.4 M , with the constraint R 1.4 = 11.75 +0.86 −0.81 km that was derived in [1] from available multi-messenger observations. This constraint cannot be fulfilled by the models based on M g = 500 MeV.
Increasing the effective gluon mass diminishes the reduction of the vector and diquark couplings at high densities. This leads to a stiffening of the quark matter EoS and increases the maximum mass of the corresponding NS sequences. As we have shown in the previous section, increasing the gluon mass lowers the NS mass where the onset of deconfinement occurs. Therefore, these masses appear to be anticorrelated.
Lowering the onset mass of quark matter for the models with M g ≥ 600 MeV also provides agreement of our approach with the constraint on R 1.4 . Thus, for all M g ≥ 600 MeV the vector and diquark couplings can be adjusted so that all the above constraints on the mass-radius relation of NSs are respected. We also would like to stress that the best agreement with the observational data is provided for low quark matter onset masses. A remarkable feature of the mass-radius relation of NSs with quark cores corresponds to existence of the so-called "special point" (SP) being a narrow region of intersection of the mass-radius curves [52]. In Refs. [52][53][54][55][56] it was reported and thoroughly studied within the CSS parameterization of the quark-matter EoS. As is seen from Fig. 4, the SP also appears within the present approach, which can not be phenomenologically caught by the CSS parameterization. Therefore, we conclude that the SP is likely to be a rather general feature of solutions of the TOV equation not limited to a given class of EoS of hybrid NSs.
The analysis of the gravitational wave signal from the inspiral phase of the binary neutron star merger gives a direct access to the dimensionless tidal deformability Λ = 2 3 k 2 C −5 expressed through the second Love number k 2 and the stellar compactness C = M/R. The measurement of the signal from the merger event GW170817 allowed to constrain the tidal deformability of a 1.4 M NS to the range Λ 1.4 = 190 +390 −120 [7]. Fig. 5 shows the dependence of Λ on M. At small values of the non-perturbative gluon mass the agreement with the observational constraint is never achieved due to late onset of quark matter. At larger values of M g it is always possible to adjust the vector and diquark couplings in order to provide a tidal deformability in the range Λ = 70 − 580. Similar to the case of the mass-radius relation, the observational data prefer the values of the vector and diquark couplings corresponding to an early onset of quark matter.
The possibility of conformal or near conformal behavior of matter in the cores of heavy NSs was recently discussed in Refs. [37,42]. We check this possibility within our model, which by construction respects the conformal limit at asymptotically high densities. As it was noticed above, the smaller the non-perturbative gluon mass the earlier the conformality is reached. Therefore, we present an analysis for M g = 600 MeV. We also consider only those parameterizations of the hybrid EoS labeled by pairs of numbers (η V , η D ), which provide consistency with the above mentioned constraints on the NS mass-radius relation and tidal deformability. Fig. 6 shows profiles of energy density and squared speed of sound for NSs of several masses obtained for the EoSs fulfilling these selection criteria. The highest NS masses correspond to the maximum ones supported by the corresponding hybrid EoSs. Both ε and c 2 S exhibit a discontinuous change at the sharp interface between quark core and hadron envelope due to the strong first order phase transition. The mixed quark-hadron phase is absent since its pressure remains constant for increasing density so that there is no pressure gradient that could balance the gravitational force which compresses the matter. We note that at M = 0.6 M for the EoS with η V = 0.330 and η D = 0.750 (left panels of Fig. 6) ε and c 2 S are continuous since quark matter does not occur in this stellar configuration. In the case of softer EoS (left panels) the energy density reached in the center of the heaviest stars is higher compared to the case of stiffer EoS (right panels). In both cases this central ε does not exceed 1200 MeV fm −3 , which is well below the range of energy densities where quark matter approaches the conformal limit.
It is important to note, that at higher M g > 600 MeV the central value of the energy density is expected to be even smaller due to a stiffening of the quark matter EOS caused by the slower melting of the vector and diquark couplings. The behavior of c 2 S supports the conclusion that quark matter inside the hybrid NSs does not reach the conformal limit. Indeed, despite decreasing towards the NS center, c 2 S does not get much smaller than 0.5 even in the cores of the heaviest NSs. Now a comment with respect to the compression modulus is in order. It quantifies the curvature of the density dependent energy per nucleon E/A = ε/n B − m N of nuclear matter, where m N is the nucleon mass. Thus, the compression modulus reads [57,58] Using the thermodynamic identity p = n 2 B ∂(ε/n B )/∂n B , we can replace the derivative ∂(ε/n B )/∂n B and arrive at K NM = 9(∂p/∂n B − 2p/n B ). Furthermore, the density derivative of the pressure in this relation can be replaced using the relation c 2 S = µ −1 B ∂p/∂n B , which is provided by the thermodynamic identity µ B = ∂ε/∂n B and definition of c 2 S . Finally, utilizing n B = (p + ε)/µ B we obtain where in the second step the pressure was expressed as p = (1/3 − δ)ε. Eq. (29) relates the compression modulus to the speed of sound and interaction measure, making K N M a quantity useful for analysing the possibility of reaching the conformal limit in NSs. At small densities, where quark matter has δ 1/3 (see Section III), the compression modulus attains a positive value K NM 9µ B c 2 S . This signals about a convex energy per baryon E/A = ε/n B − m N . At high densities, close to the conformal limit c 2 S → 1/3, δ → 0 and K NM → −3µ B /2 indicating a concave ε/n B ∝ n 1/3 B . This scaling follows from the fact that in the conformal regime ε ∝ µ 4 B and n B ∝ µ 3 B . Thus, K NM is a monotonously decreasing function of density changing from positive to negative values. Its vanishing, i.e. K NM = 0, indicates the inflection point of ε/n B , which is a precursor of the conformal regime providing K NM < 0. Fig. 7 shows profiles of the compression modulus for several NS masses obtained for the hybrid EoSs with M g = 600 MeV, which respect the observational constraints mentioned above. These profiles exhibit a discontinuous change of K NM due to the first order phase transition. It is seen from Fig. 2 that for the considered EoSs p ε at the phase transition leading to K N M 9µ B c 2 S . This explains the high values of the compression modulus of quark matter in the vicinity of the deconfinement transition since the corresponding c 2 S is almost an order of magnitude larger than in the case of nuclear matter at the saturation density. At small M, the energy density of NS matter remains small and quark part of the profiles is quite flat. In the case of heavy NSs the range of ε extends to higher values and decrease of the compression modulus toward the NS center becomes prominent. This decrease is more pronounced for softer EoSs. In other words, the smaller the NS maximum mass provided by a given hybrid EoS, the smaller values of the compression modulus are reached in its center. At the same time, observational constraints on the NS mass-radius relation require quite stiff EoS. As a result, K NM barely vanishes, if at all, even in the centers of the heaviest NSs. For example, the heaviest stellar configuration presented on the left panel of Fig. 7 yields K N M = −49 MeV, which is negligible compared to the value of 4.9 GeV on the quark matter side of the quark deconfinement transition. Thus, within the range of densities typical for NSs ε/n B remains convex or marginally reaches the inflection point meaning that quark matter remains far from the conformal limit.

V. DISCUSSION AND CONCLUSIONS
The problem of restoring the conformal limit within effective models of quark matter was given a treatment within the recently proposed relativistic density functional approach. In addition to mimicking quark confinement by a rapid growth of the quark self-energy in the confining region, the pseudo-scalar sector of this approach is equivalent to a chiral quark model with medium dependent coupling constants. In order to provide its conformal behavior at high densities we generalized the vector repulsion and diquark pairing channels to the case when the corresponding couplings decrease with the density. We also demonstrated that a quark model with density dependent vector and diquark couplings can be formulated within the relativistic density functional approach. The particular behavior of the vector and diquark couplings was motivated by an analysis of the quark repulsion energy due to non-perturbative gluon exchange in QCD in the Landau gauge. We showed that the conformal limit is asymptotically reached within the present approach at any finite value of the non-perturbative gluon mass.
The developed approach was applied in order to Maxwell-construct a family of hybrid quark-hadron EoS used to modelling NS. We found that observational data prefer the non-perturbative gluon mass exceeding a value above 500 MeV. Another general conclusion of our analysis is that color superconductivity lowers the onset density of quark matter and that such early quark deconfinement is favoured by the observational constraints on the mass-radius relation and tidal deformability of NSs. More precisely, our analysis supports the hadron-to-quark matter transition at energy densities within the range 180-500 MeV fm −3 . We also report that energy density reached in the cores of the heaviest NSs is far from the region of conformality of quark matter.
where we used definition of the quark mean-field self-energy given by Eq. (5). The fourth term in the expression for the thermodynamic potential (12) is a function of ω and q + q , entering it through the vector coupling. Thus Similarly, for the fifth term we obtain ∂G D ∂| q c iτ 2 γ 5 λ 2 q | d| q c iτ 2 γ 5 λ 2 q |.
The next step corresponds to considering the total differential of the thermodynamic potential at constant quark chemical potentials (dµ u = dµ d = 0). We note that in this case the differentials d q + q and d| q c iτ 2 γ 5 λ 2 q | do not necessarily vanish due to variation of chiral condensate, vector field and diquark pairing gap. Thus, using Eqs. (A1) -(A6) differential of the thermodynamic potential can be written as ∂G D ∂| q c iτ 2 γ 5 λ 2 q | d| q c iτ 2 γ 5 λ 2 q |.
The conditions to minimize the thermodynamic potential can be found by requiring zero values of the coefficients near the differentials of the above mentioned variables of Ω. The second line of Eq. (A7) yields the mean-field equation for the vector field, i.e. ω = −2G V q + q . With this result and q + q = f f + f , we conclude from the first line of Eq. (A7) that f + f = −∂Ω q /∂µ f . Direct differentiation of Eq. (13) with respect to µ f yields the expression for number density of a given quark flavor (18). It is seen from the third line of Eq. (A7) that the mean-field equation for chiral condensate is qq = ∂Ω q /∂m. One arrives at Eq. (19) by directly calculating the partial derivative of Ω q with respect to the current quark mass. From the fourth and fifth lines of the expression for dΩ we immediately recover the pairing gap equation ∆ = 2G D | q c iτ 2 γ 5 λ 2 q | and the diquark condensate | q c iτ 2 γ 5 λ 2 q | = −∂Ω q /∂∆. The latter can be given the form of Eq. (20) by finding the partial derivative of the thermodynamic potential of quark quasiparticles with respect to the diquark paining gap.