# Density functional approach to quark matter with confinement and color
superconductivity

Oleksii Ivanytskyi<EMAIL_ADDRESS>Institute of Theoretical
Physics, University of Wroclaw, Max Born Pl. 9, 50-204 Wroclaw, Poland  David
Blaschke<EMAIL_ADDRESS>Institute of Theoretical Physics,
University of Wroclaw, Max Born Pl. 9, 50-204 Wroclaw, Poland

###### Abstract

We present a novel relativistic density-functional approach to modeling quark
matter with a mechanism to mimic confinement. The quasiparticle treatment of
quarks provides their suppression due to a large quark selfenergy already at
the mean-field level. We demonstrate that our approach is equivalent to a
chiral quark model with medium-dependent couplings. The dynamical restoration
of the chiral symmetry is ensured by construction of the density functional.
Supplemented with the vector repulsion and diquark pairing, the model is
applied to construct a hybrid quark-hadron EoS of cold compact-star matter. We
study the connection of such a hybrid EoS with the stellar mass-radius
relation and tidal deformability. The model results are compared to various
observational constraints including the NICER radius measurement of PSR
J0740+6620 and the tidal deformability constraint from GW170817. As a striking
result we present selected parametrizations for which the hybrid star sequence
reaches masses in the range 2.5 - 2.67 $M_{\odot}$ of the lighter object in
the merger GW190814, which therefore might be interpreted as the heaviest
neutron star with to date.

###### pacs:

97.60.Jd, 26.60.Kp, 12.39.Ki

##  I Introduction

The quest for a reliable description of the deconfinement transition from
hadronic to quark matter is highly actual. A description from first principles
is possible by simulations of quantum chromodynamics (QCD) as a lattice gauge
theory, but at present it is restricted to the case of small baryochemical
potentials $\mu\lesssim 2.5~{}T$ Bazavov et al. (2019), where $T$ is the
temperature.

In order to address the modern challenges of a description of the equation of
state of matter created in the beam energy scan programs of present heavy-ion
collision (HIC) experiments at CERN and RHIC as well as those planned at FAIR
and NICA, and in the interiors of neutron stars as well as their mergers, one
needs to develop effective approaches. These should be guided by the lattice
QCD studies, neutron star phenomenology and HIC experiments. A wide class of
such models are the so-called two-phase approaches which model the hadronic
and the quark matter phases separately and then perform constructions to
describe the transitions between them.

Prominent examples, in particular for hybrid neutron star matter, are the
Maxwell construction for locally conserved charges and the Glendenning
construction Glendenning (1992) for several, globally conserved charges. Both
are based on the Gibbs principles of phase equilibrium. They can be unified by
invoking the formation of finite-size structures (pasta phases) in the
coexistence phase of quark and hadronic matter Heiselberg et al. (1993);
Yasutake et al. (2014); Maslov et al. (2019). The application of these
constructions assumes that the EoS models for the two phases are applicable in
the region of thermodynamic parameters where the Gibbs conditions for the
transition are fulfilled. Should these conditions not be fulfilled or not make
sense, this is a good signature for the inapplicability of the corresponding
EoS models in that region. In that case, interpolation constructions have
recently been suggested Masuda et al. (2013); Alvarez-Castillo et al. (2014);
Baym et al. (2018); Blaschke and Alvarez-Castillo (2020); Abgaryan et al.
(2018); Ayriyan et al. (2021); Ivanytskyi and Blaschke (2022). Such an
agnostic approach can be understood as a shortcut subsuming three steps in the
description of the hadron-to-quark matter transition: 1) the account for the
quark substructure of hadrons which leads to quark exchange interactions among
hadrons (quark Pauli blocking), 2) the account for confining interactions and
color superconductivity and 3) pasta phases in the coexistence region of quark
and hadronic matter phases.

The present work is devoted to investigating the consequences of the inclusion
of diquark interactions and thus the formation of a color superconducting
phase on the structure of the EoS and hybrid neutron star properties in a
generalization of the density functional approach Kaltenborn et al. (2017)
with a confining mean field that was motivated by the string-flip model (SFM)
Horowitz et al. (1985); Röpke et al. (1986) of quark matter. We would like to
stress, that strangeness in the density functional approach is an interesting
question, which deserves a separate study. The first steps in this direction
have been done recently Blaschke et al. (2022). However, even with the
simplifications adopted in that study, the consideration remains quite
complicated in the presence of color superconductivity. Since for a large
class of chiral quark models at zero temperature strange quarks are reported
to appear sequentially, after the light ones Klähn et al. (2007, 2013), we use
this as an argument to justify the neglect of strange quarks at the hadron-to-
quark matter transition in the present work. This assumption also keeps the
considerations sufficiently transparent.

The paper is organized as follows. A general description of the relativistic
density functional approach to quark matter within the mean field
approximation is given in the next section. Section III presents a treatment
of the pseudoscalar mode needed to determine parameters of the model with the
known mass and decay constant of the pion that is obtained within a Gaussian
approximation to the expansion beyond the mean field. The density functional
of the confining interaction potential used in this work is described in
Section IV. We discuss normalization of the model in Section V. In Section VI
we give a treatment of cold quark matter. In Section VII we construct a hybrid
quark-hadron EoS and apply it for modeling compact stars with quark matter
cores. The conclusions are given in Section VIII.

##  II Relativistic density functional approach

A general scheme to treat quark matter within the density functional formalism
under the mean-field approximation has been presented in Ref. Kaltenborn et
al. (2017). In the case of two quark flavors and three colors the effective
interaction in quark matter is accounted for by the potential $\mathcal{U}$
given in terms of quark bilinears of the form $\overline{q}~{}\hat{\Gamma}q$.
Here $q=(u,d)^{T}$ is a quark spinor and the vertex operator $\hat{\Gamma}$
acts in color, flavor and Dirac space. In the present study we attribute
scalar and pseudoscalar interaction channels to the density functional
$\mathcal{U}$. The chiral invariance of the quark interaction is ensured by
requiring $\mathcal{U}$ to depend on a chirally symmetric argument. In the two
flavor case the latter can be chosen as
$(\overline{q}q)^{2}+(\overline{q}i\gamma_{5}\vec{\tau}q)^{2}$. Clearly,
$\mathcal{U}$ can include a dependence on other quark bilinears as it was
considered in Ref. Kaltenborn et al. (2017). These complications, however, are
not necessary to maintain the proper chiral dynamics of quark matter. Due to
this we prefer to consider the simplest $\mathcal{U}$ with the above mentioned
chirally symmetric argument. Vector and diquark pairing interaction channels,
which are important for the compact star phenomenology, are incorporated to
the model within the standard framework Baym et al. (2018). Thus, the total
Lagrangian reads

$\displaystyle\mathcal{L}$ $\displaystyle=$
$\displaystyle\overline{q}(i\not{\partial}-\hat{m})q-\mathcal{U}$ (1)
$\displaystyle-$ $\displaystyle
G_{V}(\overline{q}\gamma_{\mu}q)^{2}+G_{D}(\overline{q}i\gamma_{5}\tau_{2}\lambda_{A}q^{c})(\overline{q}^{c}i\gamma_{5}\tau_{2}\lambda_{A}q),\quad$

where the diagonal matrix $\hat{m}={\rm diag}(m_{u},m_{d})$ represents current
masses of the corresponding states. Below they are labeled with the subscript
index $f=u,d$. The charge conjugate quark field is given by
$q^{c}=i\gamma_{2}\gamma_{0}\overline{q}^{T}$, while summation over the color
index $A$ in the last term of Eq. (1) and elsewhere in this paper is performed
for the antisymmetric Gell-Mann matrices $\lambda_{A}$ with $A=2,5,7$. The
strength of the interaction in vector and diquark channels is controlled by
the couplings $G_{V}$ and $G_{D}$, respectively, which are treated as free
parameters of the approach.

A Taylor expansion of the interaction potential around the expectation values
of the corresponding quark bilinears gives an efficient way of treating the
present density functional theory. The linear order of the expansion
corresponds to the selfconsistent mean-field level Kaltenborn et al. (2017),
which is applied to all interaction channels considered in this work. However,
as an additional constraint to the density functional $\mathcal{U}$ that was
not considered in Kaltenborn et al. (2017), we require the present model to
reproduce the mass and decay constant of the pion in accordance with QCD
vacuum phenomenology, thus going beyond the mean-field level in the scalar-
pseudoscalar sector. The corresponding mesonic correlations can be introduced
to the model by expanding $\mathcal{U}$ up to the second order in deviations
of quark bilinears from their mean field expectation values
$\langle\overline{q}q\rangle$ and
$\langle\overline{q}i\gamma_{5}\vec{\tau}q\rangle=0$. Below all the quantities
defined at these expectation values are denoted with the subscript index
$``MF"$. Then, the second order expanded potential energy density functional
obtains the form

$\displaystyle\mathcal{U}^{(2)}$ $\displaystyle=$
$\displaystyle\mathcal{U}_{MF}+\left(\overline{q}q-\langle\overline{q}q\rangle\right)\Sigma_{MF}$
(2) $\displaystyle-$ $\displaystyle
G_{S}(\overline{q}q-\langle\overline{q}q\rangle)^{2}-G_{PS}(\overline{q}i\gamma_{5}\vec{\tau}q)^{2},$

where the expansion coefficients are

$\displaystyle\Sigma_{MF}$ $\displaystyle=$
$\displaystyle\frac{\partial\mathcal{U}_{MF}}{\partial\langle\overline{q}q\rangle},$
(3) $\displaystyle G_{S}$ $\displaystyle=$
$\displaystyle-\frac{1}{2}\frac{\partial^{2}\mathcal{U}_{MF}}{\partial\langle\overline{q}q\rangle^{2}},$
(4) $\displaystyle G_{PS}$ $\displaystyle=$
$\displaystyle-\frac{1}{6}\frac{\partial^{2}\mathcal{U}_{MF}}{\partial\langle\overline{q}i\gamma_{5}\vec{\tau}q\rangle^{2}}.$
(5)

They represent the mean-field part of the quark self-energy $\Sigma_{MF}$,
while $G_{S}$ and $G_{PS}$ are medium dependent effective couplings in the
scalar and pseudoscalar channels, respectively. The mean field self-energy of
quarks renormalizes their mass resulting in the medium-dependent effective
mass $m^{*}_{f}=m_{f}+\Sigma_{MF}$. It is important to note that the scalar
$G_{S}$ and pseudoscalar $G_{PS}$ couplings coincide up to constant positive
factors with the second derivatives of the Lagrangian $\mathcal{L}$ with
respect to scalar $\langle\overline{q}q\rangle$ and pseudoscalar
$\langle\overline{q}i\gamma_{5}\vec{\tau}q\rangle$ condensates. Positiveness
of these derivatives provides that mean-field solution corresponds to a
minimum of the thermodynamic potential. Thus, we require $G_{S}>0$ and
$G_{PS}>0$.

The Lagrangian (1) with the density functional expanded up to the second order
terms (2) yields
$\mathcal{L}^{(2)}=\mathcal{L}_{\mathcal{U}=0}-\mathcal{U}^{(2)}$, which is
quadratic in all quark bilinears present in the model. With this Lagrangian
the partition function of the present model can be given the form of the
functional integral over the quark fields

$\displaystyle\mathcal{Z}=\int\mathcal{D}\overline{q}~{}\mathcal{D}q\exp\left[\int
dx_{E}~{}\left(\mathcal{L}^{(2)}+q^{+}\hat{\mu}q\right)\right].$ (6)

Hereafter integration over the Euclidean space-time $dx_{E}=d\tau~{}d\bf x$ is
limited to the system volume $V=\int d\bf x$ and inverse temperature
$\frac{1}{T}\equiv\beta=\int_{0}^{\beta}d\tau$. The diagonal matrix
$\hat{\mu}={\rm diag}(\mu_{u},\mu_{d})$ acting in the flavor space represents
the quark chemical potentials $\mu_{f}=\frac{\mu_{B}}{3}+Q_{f}\mu_{Q}$. They
are expressed in terms of the baryonic ($\mu_{B}$) and electric ($\mu_{Q}$)
chemical potentials as well as the baryonic (1/3) and electric ($Q_{f}$)
charges of the quarks. The partition function provides direct access to the
thermodynamic potential $\Omega=-(T/V)\ln\mathcal{Z}$ from which all equations
of state can be derived.

To proceed further we bosonize $\mathcal{L}^{(2)}$ by means of the Hubbard-
Stratonovich transformation in the spirit of Ref. Blaschke et al. (2014). This
introduces collective scalar $\sigma$, pseudoscalar $\vec{\pi}$, vector
$\omega_{\mu}$ and complex scalar diquark $\Delta_{A}$ fields, coupled to
$\overline{q}q-\langle\overline{q}q\rangle$,
$\overline{q}i\gamma_{5}\vec{\tau}q$, $\overline{q}\gamma_{\mu}q$ and
$\overline{q}i\gamma_{5}\tau_{2}\lambda_{A}q$, respectively. Then the
partition function becomes

$\displaystyle\mathcal{Z}$ $\displaystyle=$
$\displaystyle\int\mathcal{D}\sigma~{}\mathcal{D}\vec{\pi}~{}\mathcal{D}\omega_{\mu}~{}\mathcal{D}\Delta_{A}~{}\mathcal{D}\Delta_{A}^{*}~{}\mathcal{D}\overline{q}~{}\mathcal{D}q$
(7) $\displaystyle\times$ $\displaystyle\exp\left[\int
dx_{E}\left(\mathcal{L}^{bos}+q^{+}\hat{\mu}q\right)\right].$

The bosonized Lagrangian $\mathcal{L}^{bos}$ can be expressed through the
Nambu-Gorkov bispinor quark fields
$\mathcal{Q}^{T}=\frac{1}{\sqrt{2}}(q~{}q^{c})$ as

$\displaystyle\mathcal{L}^{bos}+q^{+}\hat{\mu}q$ $\displaystyle=$
$\displaystyle\overline{\mathcal{Q}}\mathcal{S}^{-1}\mathcal{Q}-\mathcal{U}_{MF}+\langle\overline{q}q\rangle\Sigma_{MF}+\sigma\langle\overline{q}q\rangle$
(8) $\displaystyle-$
$\displaystyle\frac{\sigma^{2}}{4G_{S}}-\frac{\vec{\pi}^{2}}{4G_{PS}}+\frac{\omega_{\mu}\omega^{\mu}}{4G_{V}}-\frac{\Delta_{A}^{*}\Delta_{A}}{4G_{D}}.$

Here the inverse Nambu-Gorkov quark propagator is

$\displaystyle\mathcal{S}^{-1}=\left(\begin{array}[]{l}S^{-1}_{+}-\sigma-i\gamma_{5}\vec{\tau}\cdot\vec{\pi}\hskip
28.45274pti\Delta_{A}\gamma_{5}\tau_{2}\lambda_{A}\\\ \hskip
14.22636pti\Delta_{A}^{*}\gamma_{5}\tau_{2}\lambda_{A}\hskip
28.45274ptS^{-1}_{-}-\sigma-i\gamma_{5}\vec{\tau}^{T}\cdot\vec{\pi}\end{array}\right),$
(11)

with
$S^{-1}_{\pm}=i\not{\partial}\pm\not{\omega}-m^{*}\pm\gamma_{0}\hat{\mu}$. The
quark fields in Eq. (7) can be integrated out since $\mathcal{L}^{bos}$ is
quadratic with respect to them. This yields the partition function

$\displaystyle\mathcal{Z}$ $\displaystyle=$
$\displaystyle\int\mathcal{D}\sigma~{}\mathcal{D}\vec{\pi}~{}\mathcal{D}\omega_{\mu}~{}\mathcal{D}\Delta_{A}~{}\mathcal{D}\Delta_{A}^{*}$
(12) $\displaystyle\exp\left[\frac{1}{2}{\rm
Tr}\ln(\beta\mathcal{S}^{-1})-\int
dx_{E}\Biggl{(}\mathcal{U}_{MF}-\langle\overline{q}q\rangle\Sigma_{MF}\right.$
$\displaystyle\left.-\sigma\langle\overline{q}q\rangle+\frac{\sigma^{2}}{4G_{S}}+\frac{\vec{\pi}^{2}}{4G_{PS}}-\frac{\omega_{\mu}\omega^{\mu}}{4G_{V}}+\frac{\Delta_{A}^{*}\Delta_{A}}{4G_{D}}\Biggl{)}\right].$

From now on we prefer to work in the momentum representation rather than in
the coordinate one. Thus, for the quark propagator in Eq. (11) we use the
Fourier transformed bosonic fields and $S^{-1}_{\pm}=\not{k}-\hat{m}^{*}$ with
$k_{0}=iz_{n}\pm\hat{\mu}^{*}$ and $z_{n}=(2n+1)\pi T$ being a fermionic
Matsubara frequency. The trace in Eq. (12) is performed in the Nambu-Gorkov,
Dirac, flavor, and color spaces, as well as over three-momentum and Matsubara
indices. We note that the factor $\frac{1}{2}$ compensates the artificial
doubling of quark degrees of freedom in the Nambu-Gorkov formalism.

This work is focused on the mean field treatment of the present model. In this
case the functional integration over the bosonic fields is dropped, while the
fields themselves are replaced by their values at the extremum of the
Euclidean action. These expectation values can be deduced from the gap
equations obtained by averaging the corresponding Euler-Lagrange equations,

$\displaystyle\sigma$ $\displaystyle=$
$\displaystyle-2G_{S}(\overline{q}q-\langle\overline{q}q\rangle),$ (13)
$\displaystyle\vec{\pi}$ $\displaystyle=$
$\displaystyle-2G_{PS}\overline{q}i\gamma_{5}\vec{\tau}q,$ (14)
$\displaystyle\omega_{\mu}$ $\displaystyle=$
$\displaystyle-2G_{V}\overline{q}\gamma_{\mu}q,$ (15)
$\displaystyle\Delta_{A}$ $\displaystyle=$ $\displaystyle
2G_{D}\overline{q}^{c}i\tau_{2}\gamma_{5}\lambda_{A}q.$ (16)

We note that the equation for $\Delta_{A}^{*}$ is equivalent to the hermitean
conjugate of Eq. (16). It follows from the equations for (pseudo)scalar fields
that $\langle\sigma\rangle=\langle\vec{\pi}\rangle=0$ under the mean field
approximation. Therefore, the contribution of $\sigma$ and $\vec{\pi}$ to the
thermodynamic properties of the present model is neglected at this level.
Below, however, we study the pseudoscalar mode within the Gaussian
approximation in order to normalize our approach to the vacuum phenomenology
of QCD. The expectation values of the vector and diquark fields can be
obtained by averaging Eqs. (15) and (16), respectively. By a proper Lorentz
transform, the vector field attains the form
$\langle\omega_{\mu}\rangle=g_{\mu 0}\omega$. The constant
$\omega=-2G_{V}\langle q^{+}q\rangle$ is absorbed to the effective chemical
potential of quarks $\mu^{*}_{f}=\mu_{f}+\omega$. Furthermore, we perform a
global color rotation leaving $\Delta_{2}$ the only diquark field with non
vanishing expectation value. Only its modulus
$\Delta=2G_{D}|\langle\overline{q}^{c}i\tau_{2}\gamma_{5}\lambda_{2}q\rangle|$
appears in the expressions for thermodynamic quantities. This non-vanishing
diquark field is conjugated to Gell-Mann matrix $\lambda_{2}$ pairing only red
and green quark states, while leaving the blue one unpaired. We label these
states with the subscript index $c=(r,g,b)$. The corresponding single particle
energies shifted by the effective chemical potential $\mu_{f}^{*}$ can be
found from the condition $\det(\mathcal{S}^{-1}_{MF})=0$ as

$\displaystyle\epsilon^{\pm}_{{\bf k}fc}={\rm sgn}(\epsilon_{{\bf
k}f}\mp\mu_{f}^{*})\sqrt{(\epsilon_{{\bf
k}f}\mp\mu_{f}^{*})^{2}+\Delta^{2}_{c}},$ (17)

where the subscript index $\bf k$ labels quark momentum states,
$\epsilon_{{\bf k}f}=\sqrt{{\bf k}^{2}+m_{f}^{*2}}$ and
$\Delta_{r}=\Delta_{g}=\Delta$, $\Delta_{b}=0$ are introduced for unifying the
notations. Hereafter, the superscript indices $``+"$ and $``-"$ correspond to
quarks and antiquarks, respectively. It follows from this expression that
finite $\Delta$ generates a gap in the dispersion relation of red and green
quarks at $\mu_{f}^{*}>m_{f}^{*}$ and antiquarks of the same colors at
$\mu_{f}^{*}<-m_{f}^{*}$.

With the above notations, the mean-field thermodynamic potential can be
written as

$\displaystyle\Omega=\Omega_{q}+\mathcal{U}_{MF}-\langle\overline{q}q\rangle\Sigma_{MF}-\frac{\omega^{2}}{4G_{V}}+\frac{\Delta^{2}}{4G_{D}}.$
(18)

Its quark part

$\displaystyle\Omega_{q}$ $\displaystyle=$ $\displaystyle-\frac{T}{2V}{\rm
Tr}\ln(\beta\mathcal{S}^{-1}_{MF})$ (19) $\displaystyle=$
$\displaystyle-2\sum_{f,c,a=\pm}\int\frac{d{\bf
k}}{(2\pi)^{3}}\left[\frac{g_{\bf k}}{2}\epsilon_{{\bf
k}fc}^{a}-T\ln\left(1-f^{a}_{{\bf k}fc}\right)\right]$

includes the single particle distribution function $f_{{\bf
k}fc}^{\pm}=\left[e^{\beta\epsilon_{{\bf k}fc}^{\pm}}+1\right]^{-1}$. The zero
point terms in Eq. (19) are regularized by the formfactor

$\displaystyle g_{\bf k}=\exp\left[-\left(\frac{\bf
k^{2}}{\Lambda^{2}}\right)^{\xi}\right]$ (20)

with constant $\xi>{1}/{2}$ and $\Lambda$ defining a momentum scale. This
formfactor has an inflection point at $|{\bf
k}|=\Lambda\left(1-\frac{1}{2\xi}\right)^{1/{2\xi}}$. In this work we use its
Gaussian form with $\xi=1$. Regularization with a sharp momentum cutoff
corresponds to $\xi\rightarrow\infty$ yielding $g_{\bf
k}\rightarrow\theta(\Lambda-|{\bf k}|)$. However, we omit this sharp
regularization scheme since it leads to unphysical irregularities in the
behavior of various thermodynamic quantities of the 2SC phase. In order to
show this we consider the number density $\langle
f^{+}f\rangle=-{\partial\Omega_{q}}/{\partial\mu_{f}}$ for a given quark
flavor $f$. Differentiation of ${\rm sgn}(\epsilon_{{\bf k}f}\mp\mu_{f}^{*})$
in the dispersion relation of (anti)quarks produces the Dirac delta-function
$\delta(\epsilon_{{\bf k}f}\mp\mu_{f}^{*})$. The part of $\langle
f^{+}f\rangle$, which include this delta-function, is

$\displaystyle\delta\langle f^{+}f\rangle$ $\displaystyle=$ $\displaystyle
2\sum_{c,a=\pm}a\int\frac{d{\bf k}}{(2\pi)^{3}}\left(f_{{\bf
k}fc}^{a}-\frac{g_{\bf k}}{2}\right)2\Delta_{c}\delta(\epsilon_{{\bf
k}fb}^{a})$ (21) $\displaystyle=$
$\displaystyle\frac{2\mu_{f}^{*}k_{f}\Delta}{\pi^{2}}(1-g_{\bf k})_{|{\bf
k}|=k_{f}},$

where
$k_{f}=\theta\left(\mu_{f}^{*2}-m_{f}^{*2}\right)\sqrt{\mu_{f}^{*2}-m_{f}^{*2}}$
is the Fermi momentum of a given quark flavor. From this expression it is
clear that any discontinuity of the formfactor $g_{\bf k}$ translates to
behavior of quark number density if $\Delta\neq 0$. The same conclusion holds
for any thermodynamic quantity of the 2SC phase. This explains the choice of
the momentum formfactor $g_{\bf k}$.

The mean-field equation for the zeroth component of vector field is

$\displaystyle\omega$ $\displaystyle=$ $\displaystyle
2G_{V}\sum_{f}\frac{\partial\Omega_{q}}{\partial\mu_{f}}$ $\displaystyle=$
$\displaystyle 4G_{V}\sum_{f,c,a=\pm}a\int\frac{d{\bf
k}}{(2\pi)^{3}}\left(\frac{g_{\bf k}}{2}-f_{{\bf
k}fc}^{a}\right)\left(2\Delta_{c}\delta(\epsilon_{{\bf
k}fb}^{a})+\frac{\epsilon_{{\bf k}fb}^{a}}{\epsilon_{{\bf k}fc}^{a}}\right).$

The pairing gap equation can be derived in a similar way

$\displaystyle\Delta$ $\displaystyle=$
$\displaystyle-2G_{D}\sum_{f}\frac{\partial\Omega_{q}}{\partial\Delta}$ (23)
$\displaystyle=$ $\displaystyle 4G_{D}\sum_{f,c,a=\pm}\int\frac{d{\bf
k}}{(2\pi)^{3}}\left(\frac{g_{\bf k}}{2}-f_{{\bf
k}fc}^{a}\right)\frac{\Delta_{c}}{\epsilon_{{\bf k}fc}^{a}}.$

This equation has two solutions. The trivial one, $\Delta=0$, corresponds to
the normal phase. Under suitable conditions, Eq. (23) also has a nontrivial
solution $\Delta\neq 0$, signaling the formation of a color superconducting
phase of quark matter. Due to the presence of an attractive quark-quark
interaction, according to the Cooper theorem Cooper (1956), the occurrence of
such a phase is inevitable at sufficiently low temperatures and high quark
number densities. Finally, the equation for the chiral condensate reads

$\displaystyle\langle\overline{q}q\rangle$ $\displaystyle=$
$\displaystyle\sum_{f}\frac{\partial\Omega_{q}}{\partial m_{f}}$
$\displaystyle=$ $\displaystyle 2\sum_{f,c,a=\pm}\int\frac{d{\bf
k}}{(2\pi)^{3}}\left(f_{{\bf k}fc}^{a}-\frac{g_{\bf
k}}{2}\right)\left(2\Delta_{c}\delta(\epsilon_{{\bf
k}fb}^{a})+\frac{\epsilon_{{\bf k}fb}^{a}}{\epsilon_{{\bf
k}fc}^{a}}\right)\frac{m^{*}}{\epsilon_{{\bf k}f}}.$

Eqs. (II)-(II) provide the stationarity of the thermodynamic potential. Having
these conditions solved, the pressure, the number density of a given quark
flavor, the entropy density and the energy density can be found using standard
thermodynamic identities such as $p=-\Omega+\Omega_{0}$, $n_{f}={\partial
p}/{\partial\mu_{f}}$, $s={\partial p}/{\partial T}$,
$\varepsilon=\sum_{f}\mu_{f}n_{f}+Ts-p$, respectively. Hereafter, the
subscript $``0"$ denotes the quantities defined in the vacuum, i.e. at
$\hat{\mu}=T=0$.

##  III Scalar and pseudoscalar modes

We now proceed to the analysis of the (pseudo)scalar fluctuation modes
$\Sigma=\Sigma_{\sigma}+\Sigma_{\pi}$ with

$\displaystyle\Sigma_{\pi}=-\left(\begin{array}[]{l}i\gamma_{5}\vec{\tau}\vec{\pi}\hskip
17.07182pt0\\\ \hskip 8.5359pt0\hskip
14.22636pti\gamma_{5}\vec{\tau}^{T}\vec{\pi}\end{array}\right),\quad\Sigma_{\sigma}=-\left(\begin{array}[]{l}\sigma\hskip
5.69046pt0\\\ 0\hskip 5.69046pt\sigma\end{array}\right),$ (29)

that appear beyond the mean field approximation as dressing of the inverse
Nambu-Gorkov quark propagator $\mathcal{S}^{-1}=\mathcal{S}^{-1}_{MF}+\Sigma$.
Within the Gaussian approximation, $\Sigma$ is assumed to be small compared to
$\mathcal{S}^{-1}_{MF}$. Thus, the logarithm of the inverse dressed quark
propagator can be expanded as

$\displaystyle\ln\left(\beta\mathcal{S}^{-1}\right)$ $\displaystyle=$
$\displaystyle\ln\left(\beta\mathcal{S}^{-1}_{MF}\right)+\ln\left(1+\mathcal{S}_{MF}\Sigma\right)$
(30) $\displaystyle=$
$\displaystyle\ln\left(\beta\mathcal{S}^{-1}_{MF}\right)+\mathcal{S}_{MF}\Sigma-\frac{\left(\mathcal{S}_{MF}\Sigma\right)^{2}}{2}+\mathcal{O}\left(\Sigma^{3}\right).\quad$

The trace of the pion contribution in the first order term in this expansion
vanishes due to the tracelessness of the flavor Pauli matrices. The trace of
$\mathcal{S}_{MF}\Sigma_{\sigma}$ exactly cancels the term $\int
dx_{E}~{}\sigma\langle\overline{q}q\rangle$ in Eq. (12) and, therefore, both
terms vanish. The pion-sigma cross contribution
$\propto\Sigma_{\sigma}\Sigma_{\pi}$ in the second order term also vanishes
due to traceless character of the Pauli matrices. The remaining second order
terms $\sim\Sigma_{\sigma}^{2}$ and $\sim\Sigma_{\pi}^{2}$ are quadratic in
the quark propagator and represent one-loop polarization operators of scalar
and pseudoscalar modes at four-momentum $p$,

$\displaystyle\Pi_{\pi}$ $\displaystyle=$ $\displaystyle-\frac{1}{2\beta
V}{\rm
Tr}(i\gamma_{5}\mathcal{S}_{MF})^{2}=-\frac{\langle\overline{q}q\rangle}{m^{*}}+p^{2}\mathcal{I}(p^{2}),$
(31) $\displaystyle\Pi_{\sigma}$ $\displaystyle=$
$\displaystyle-\frac{1}{2\beta V}{\rm
Tr}(\mathcal{S}_{MF})^{2}=-\frac{\langle\overline{q}q\rangle}{m^{*}}+(p^{2}-4m^{*2})\mathcal{I}(p^{2}).\quad$
(32)

Here the integral $\mathcal{I}(p^{2})$ has been introduced as

$\displaystyle\mathcal{I}(p^{2})=2\sum_{f,c,a=\pm}\int\frac{d{\bf
k}}{(2\pi)^{3}}\frac{g_{\bf k}-2f_{{\bf k}fc}^{a}}{2\epsilon_{{}_{\bf
k}f}(4\epsilon_{{\bf k}f}^{2}-p^{2})}.$ (33)

For the sake of simplicity we evaluated $\Pi_{\sigma}$ and $\Pi_{\pi}$ at
$\Delta=0$, which is sufficient to define the vacuum properties of the
(pseudo)scalar modes. The details of the derivation can be found in Refs.
Klevansky (1992); Blaschke et al. (2014); Hüfner et al. (1994); Zhuang et al.
(1994); Buballa (2005). These polarization operators enter the inverse of the
Fourier transformed mesonic propagators as

$\displaystyle D_{\sigma}^{-1}=\frac{1}{2G_{S}}-\Pi_{\sigma},\quad
D_{\pi}^{-1}=\frac{1}{2G_{PS}}-\Pi_{\pi}.$ (34)

These two-dimensional matrices act in the momentum-Matsubara space. The poles
of $D_{\sigma}$ and $D_{\pi}$ at ${\bf p}=0$ define mesonic masses being
solutions of

$\displaystyle M_{\sigma}^{2}$ $\displaystyle=$
$\displaystyle\left(\frac{1}{2G_{S}}+\frac{\langle\overline{q}q\rangle}{m^{*}}\right)\mathcal{I}^{-1}(M_{\sigma}^{2})+4m^{*2},$
(35) $\displaystyle M_{\pi}^{2}$ $\displaystyle=$
$\displaystyle\left(\frac{1}{2G_{PS}}+\frac{\langle\overline{q}q\rangle}{m^{*}}\right)\mathcal{I}^{-1}(M_{\pi}^{2}).$
(36)

With the effective quark mass expressed as
$m^{*}=m-2G_{PS}\langle\overline{q}q\rangle$, the equation for $M_{\pi}^{2}$
can be given a form of the Gell-Mann-Oakes-Renner relation Gell-Mann et al.
(1968) $M_{\pi}^{2}F_{\pi}^{2}=-m\langle\overline{q}q\rangle$ with the pion
decay constant

$\displaystyle F_{\pi}^{2}=m^{*}(m^{*}-m)\mathcal{I}(M_{\pi}^{2}).$ (37)

Eqs. (35) - (37) give a direct access to the vacuum phenomenology of QCD.

##  IV Density functional

The confining aspect of the quark interaction is modeled within our approach
by a large and positive scalar quark self-energy $\Sigma_{MF}$ in the
confining region. The SFM argument implies that the quark mass grows
proportional to the mean interquark separation or, equivalently,
$\Sigma\propto D_{0}(q^{+}q)^{-{1}/{3}}$ Horowitz et al. (1985); Röpke et al.
(1986), with $D_{0}$ being proportional to the string tension. It is worth
mentioning, that an alternative prescription for the low density behavior of
the effective quark mass implies $m^{*}\propto\langle q^{+}q\rangle^{-1}$
Fowler et al. (1981). For the reader’s convenience, we demonstrate in the
Appendix A how this result can be obtained under some simplifications when
compared to our approach. According to Eq. (3), the SFM behavior of $\Sigma$
corresponds to $\mathcal{U}\propto(q^{+}q)^{\frac{2}{3}}$. In the confining
region of low temperature and chemical potential, the quark number density and
the chiral condensate are connected by the approximate relation
$\langle\overline{q}q\rangle\simeq\langle\overline{q}q\rangle_{0}+\langle
q^{+}q\rangle$. This motivates
$\mathcal{U}\propto(\langle\overline{q}q\rangle_{0}-\overline{q}q)^{\frac{2}{3}}$.
At the same time, the chiral invariance of the Lagrangian (1) is ensured
within our approach by the appropriate choice of the argument
$(\overline{q}q)^{2}+(\overline{q}i\gamma_{5}\vec{\tau}q)^{2}$ in the ansatz
for the density functional $\mathcal{U}$. In order to provide such property of
the chirally symmetric interaction mimicking confinement by a significant
growth of the quark self-energy we consider the parameterization

$\displaystyle\mathcal{U}=D_{0}\left[(1+\alpha)\langle\overline{q}q\rangle_{0}^{2}-(\overline{q}q)^{2}-(\overline{q}i\gamma_{5}\vec{\tau}q)^{2}\right]^{\varkappa},$
(38)

where $D_{0}$ is the coupling strength and $\alpha$ is a constant parameter
regulating the value of the quark mass in the vacuum. Setting
$\varkappa={1}/{3}$ is in accordance to the SFM prescription adopted in
Kaltenborn et al. (2017). At zero baryonic density, Eq. (38) reproduces the
parameterization used in Ref. Kaltenborn et al. (2017) if the no-sea
approximation is applied which neglects the quark vacuum term in the
thermodynamic potential and the pseudoscalar interaction channel is neglected.
Indeed, in this approximation chiral condensate coincides up to the sign with
the quark scalar density, i.e. $n_{s}=-\langle\overline{q}q\rangle$ and
$\langle\overline{q}q\rangle_{0}=0$. This makes Eq. (38) of the present work
equivalent to Eq. (23) of Ref. Kaltenborn et al. (2017). The role of $\alpha$
becomes clear from considering effective quark mass in the vacuum when
$\langle\overline{q}q\rangle=\langle\overline{q}q\rangle_{0}$ and
$\langle\overline{q}i\gamma_{5}\vec{\tau}q\rangle=0$. With the definition
$\hat{m}^{*}=\hat{m}+\Sigma_{MF}$ and the mean field self energy of quarks
given by Eq. (3) we obtain

$\displaystyle\hat{m}^{*}_{0}=\hat{m}+\frac{2}{3}D_{0}\alpha^{-2/3}\langle\overline{q}q\rangle_{0}^{-1/3},$
(39)

At high $\alpha$ effective quark mass coincides with the current one, which
means that quark interaction in the (pseudo)scalar channels is suppressed. On
the other hand, at $\alpha=0$ vacuum effective mass diverges making any
excitations of quark degrees of freedom energetically disfavored. At small
positive values of this parameter $\hat{m}^{*}_{0}$ remains finite but large
leading to sufficient suppression of quarks. This simple physical mechanism
provides an effective phenomenological confinement within the present model.
In order to illustrate the role of $\alpha$ we show in Fig. 1 the dependence
of the effective quark mass on the baryon density found for cold symmetric
quark matter in the absence of diquark pairing. This $m^{*}$ decreases with
growing $n_{B}$ and vanishes asymptotically. Larger $\alpha$ parameters
correspond to a faster decrease of $m^{*}$ and to its smaller vacuum value,
which remains finite at all $\alpha>0$. At $\alpha=0$, the effective quark
mass in the vacuum diverges similarly to the case of the SFM Kaltenborn et al.
(2017) shown for the sake of comparison.

Eq. (38) provides $\hat{m}^{*}=\hat{m}-2G_{PS}\langle\overline{q}q\rangle$
with $G_{PS}$ defined by Eq. (5). This relation is similar to the one of the
NJL model Klevansky (1992). For $\varkappa=1$ the two models coincide. The
insignificant constant term
$-D_{0}(1+\alpha)\langle\overline{q}q\rangle_{0}^{2}$ in the expression for
the Lagrangian can be eliminated by $\alpha=-1$. In our approach, however,
effective couplings are medium dependent. As is seen from Fig. 2, they
decrease with temperature and chemical potential, which is in line with the
weakening of the strong interaction due to the running of the QCD coupling
Gross and Wilczek (1973); Politzer (1973); Deur et al. (2016). It is worth
mentioning, that a large value of $G_{PS}$ is another reason for high quark
masses in the confining region.

Furthermore, contrary to the NJL model, the couplings in the scalar and
pseudoscalar channels of the present model do not coincide in the general
case, i.e. $G_{S}\neq G_{PS}$. This signals the violation of the chiral
symmetry on the level of the Lagrangian $\mathcal{L}^{(2)}$ and is a direct
consequence of expanding $\mathcal{U}$ around the mean-field solution, which
is known to break the chiral symmetry. As can be seen from Fig. 2, at certain
values of $T$ and/or $\mu_{B}$, however, these couplings converge to the same
constant value $G_{\infty}$ as a consequence of the dynamical restoration of
the chiral symmetry. The value for $G_{\infty}$ is independent of the coupling
constants in the vector and diquark channels. For the model parameters from
Table 1 the model yields $G_{\infty}=6.6~{}{\rm GeV}^{-2}$. It is remarkable
that this asymptotic value of scalar and pseudoscalar couplings corresponds to
$G_{\infty}\Lambda^{2}=2.16$, which is very close to the value
$G_{NJL}\Lambda^{2}_{NJL}=2.14$ obtained in the NJL model Ratti et al. (2006).

Figure 1: Effective quark mass $m^{*}$ as functions of baruonic density
$n_{B}$ at different values of the parameter $\alpha$. Calculations are
performed for cold symmetric quark matter and $G_{V}=G_{D}=0$ while the rest
of the model parameters except $\alpha$ have the values from Table 1. Blue
long dotted line is obtained within the SFM with $\alpha_{SFM}=0.39~{}fm^{-3}$
(see Ref. Kaltenborn et al. (2017) for details).

##  V Parametrization of the model

$m$ [MeV] |  $\Lambda$ [MeV] | $\alpha$ | $D_{0}\Lambda^{-2}$
---|---|---|---
4.2 | 573 | 1.43 | 1.39
$M_{\pi}$ [MeV] |  $F_{\pi}$ [MeV] |  $M_{\sigma}$ [MeV] |  $\langle\overline{l}l\rangle_{0}^{1/3}$ [MeV]
140 | 92 | 980 | -267
Table 1: Parameters of the model used and resulting physical quantities.

The model has six independent parameters: the current quark mass $m$ which is
degenerate for the two flavors, the parameters of the density functional,
$\alpha$ and $D_{0}$, the momentum scale parameter $\Lambda$, the vector
coupling $G_{V}$ and the diquark coupling $G_{D}$. The first four of them
define the vacuum state of the model while the coupling constants $G_{V}$ and
$G_{D}$ are free parameters determining its high-density behavior. The
consideration of the present section is limited to the vacuum state. Most
important for the parametrization of the model are the properties of the
pseudoscalar meson, the pion. Namely its mass $M_{\pi}$ and decay constant
$F_{\pi}$, with their well-known values given in Table 1. The vacuum value of
the chiral condensate per flavor
$\langle\overline{l}l\rangle_{0}={\langle\overline{q}q\rangle_{0}}/2$ is
another important quantity. Predictions for $\langle\overline{q}q\rangle_{0}$
from chiral quark models of NJL type with a soft cutoff Grigorian (2007) like
in Eq. (II) are in tension with the result obtained from QCD sum rules at the
renormalization scale of 1 GeV, which is Jamin (2002)
$|\langle\bar{l}l\rangle_{0}^{1/3}|\big{|}_{\rm 1~{}GeV}=242\pm 15$ MeV.
Therefore, we allow a deviation of this quantity from
$\langle\overline{l}l\rangle_{0}\big{|}_{\rm 1~{}GeV}$. The scalar meson mass
$M_{\sigma}$ is also considered. Despite possibly being the lightest scalar
meson $f_{0}(500)$ has rather large width 500 - 1000 MeV Tanabashi et al.
(2018), which makes its relevance to the vacuum phenomenology quite doubtful.
Therefore, for our analysis we use the $f_{0}(980)$ meson with a width
$\lesssim$ 100 MeV Tanabashi et al. (2018). Its mass is varied within the
range from 880 MeV to 1080 MeV. We use the following strategy. For a given
value of the effective quark mass $m^{*}$ in the vacuum we fit $m$ and
$\Lambda$ to $M_{\pi}$ and $F_{\pi}$. This allows us to find
$\langle\overline{l}l\rangle_{0}$ as well as the pseudoscalar coupling
$G_{PS}=(m-m^{*})/(2\langle\overline{q}q\rangle_{0})$. Having the above
parameters fitted and varying $M_{\sigma}$ within the mentioned interval we
can calculate the scalar coupling $G_{S}$ from Eq. (35). Fig. 3 shows $G_{PS}$
and $G_{S}$ as functions of $m^{*}$. The calculations are performed for
different $m^{*}$ which provide $G_{S}>G_{PS}$ and thus represent a range of
physical values of the model parameters. Indeed, we can show by direct
calculations that in the vacuum

Figure 2: Scaled effective scalar $G_{S}\Lambda^{2}$ (thin solid curves) and
pseudoscalar $G_{PS}\Lambda^{2}$ (thick solid curve) couplings as functions of
temperature $T$ at $\mu_{B}=0$ (upper panel) and baryonic chemical potential
$\mu_{B}$ at $T=0$ (lower panel). Dotted curves on the lower panel demonstrate
unstable parts removed by applying the Maxwell construction. Dashed line
represents the NJL value $G_{NJL}\Lambda_{NJL}^{2}=2.14$ from Ref. Ratti et
al. (2006). Calculations are performed for symmetric quark matter and
$G_{V}=G_{D}=0$ while the rest of the model parameters have the values from
Table 1.

$\displaystyle
G_{S}-G_{PS}=\frac{2D_{0}}{9(1+\alpha)^{5/3}\langle\overline{q}q\rangle_{0}^{4/3}}>0.$
(40)

This quantity vanishes only for $\alpha\rightarrow\infty$ since finite scalar
and pseudoscalar couplings require $D_{0}\neq 0$. Thus, $\alpha$ remains
finite only when the effective quark mass $m^{*}$ exceeds some value, which
mildly depends on the mass of the scalar mode. It is seen from Fig. 3 that
this value grows with $M_{\sigma}$.

For a given set of the model parameters the pseudocritical temperature
$T_{PC}$ can be found according to its definition as the peak position of the
chiral susceptibility, which is defined as the temperature derivative of the
chiral condensate at vanishing chemical potential. This quantity depends on
both effective quark mass $m^{*}$ and mass of scalar meson $M_{\sigma}$ in the
vacuum. Note that $T_{PC}$ calculated in such a way is just a preliminary
result since it is limited to the approximation that is employed when
calculating the chiral condensate. When the mean-field approximation is used
then it does not account for the hadronic correlations of quarks. Such
correlations are known to contribute to the melting of the chiral condensate
in the pseudocritical region Hüfner et al. (1994); Zhuang et al. (1994) and
lead to a lowering of the pseudocritical temperature. Therefore, we require
our $T_{PC}$ to exceed the lattice QCD result $T_{PC}^{\rm lQCD}=156.5\pm 1.5$
MeV Bazavov et al. (2019). Fig. 4 shows
$|\langle\overline{l}l\rangle_{0}|^{1/3}$ and $T_{PC}$ as functions of
$m^{*}$. At large effective quark mass the chiral condensate is too high at
any value of the scalar meson mass. Therefore, we require $m^{*}\lesssim 800$
MeV. At the same time we prefer to have $m^{*}$ high enough to provide an
efficient suppression of quarks in the confinement region. This excludes low
$M_{\sigma}$, since in this case either $T_{PC}<T_{PC}^{\rm lQCD}$ or $m^{*}$
is not high enough. On the other hand, for $M_{\sigma}=1080$ MeV the
pseudocritical temperature $T_{PC}$ of the model exceeds the one from lattice
QCD $T_{PC}^{\rm lQCD}$ by more than 10%. Due to this we use the central value
of $M_{\sigma}=980$ MeV. At $m^{*}=718$ MeV this value of the scalar meson
mass leads to a reasonable value of $T_{PC}=163$ MeV and
$|\langle\overline{l}l\rangle_{0}|^{1/3}=267$ MeV compatible with results of
other chiral models Grigorian (2007). The corresponding values of the model
parameters are summarized in Table 1. It is worth mentioning that a relatively
large value of the vacuum effective quark mass in our approach is quite in
line with the results of an analysis of the data on the lattice QCD
thermodynamics performed with the quark quasiparticle model yielding
$m^{*}_{0}=610-950$ MeV Plumari et al. (2011) or even $m_{0}^{*}>1$ GeV Bozek
et al. (1998).

Figure 3: Dimensionless scalar ($G_{S}\Lambda^{2}$, thin curves) and
pseudoscalar ($G_{PS}\Lambda^{2}$, thick curve) couplings as functions of the
effective quark mass $m^{*}$ calculated for the vacuums state. Figure 4:
Chiral condensate per flavor $|\langle\overline{l}l\rangle_{0}|^{1/3}$ (thick
curve) and pseudocritical temperature $T_{PC}$ (thin curves) as functions of
the the effective quark mass $m^{*}$ in the vacuum. Filled circles represent
the low mass edges of the intervals where parameters of the model attain
physical values (see description in the text). The orange shaded area covers
the range of $T_{PC}$ generated by the considered vacuum values of $m^{*}$ and
$M_{\sigma}$. Dashed horizontal lines and grey shaded areas represent
$|\langle\overline{l}l\rangle_{0}^{\rm 1~{}GeV}|^{1/3}$ and $T_{PC}^{\rm
lQCD}$ along with the corresponding error bars obtained with QCD sum rules and
lattice QCD, respectively.

First principle information about vector and diquark couplings is absent at
the moment. The argumentation based on a Fierz transformation of the (massive)
vector boson exchange interaction model which yields $G_{V}=0.5~{}G_{S}$ and
$G_{D}=0.75~{}G_{S}$ Buballa (2005) is not directly applicable to the present
model since its (pseudo)scalar couplings are not constant. Therefore, $G_{V}$
and $G_{D}$ are treated as free parameters. We use the vacuum value of the
scalar coupling $G_{S0}=18.1~{}{\rm GeV}^{-2}$ in order to parameterize them
by $\eta_{V}\equiv\frac{G_{V}}{G_{S0}}$ and
$\eta_{D}\equiv\frac{G_{D}}{G_{S0}}$. The model parameterizations considered
below are labeled with a pair of numbers corresponding to scaled couplings
$\eta_{V}$ and $\eta_{D}$. For example, $(0.2,1.5)$ stands for a quark matter
EoS obtained for $\eta_{V}=0.2$, $\eta_{D}=1.5$ and other parameters from
Table 1.

##  VI Cold quark matter

From now on we consider the case of zero temperature, which is of interest for
the astrophysical applications. The first question we address is related to
onset of the 2SC phase. It occurs when two solutions of Eq. (23) coincide
provided by

$\displaystyle\frac{\partial^{2}\Omega}{\partial\Delta^{2}}\biggl{|}_{\Delta=0}=0.$
(41)

Solving this equation with respect to baryonic chemical potential yields its
critical value $\mu_{B}^{c}$ corresponding to the onset of color
superconductivity. Fig. 5 shows it as a function of the scaled diquark
coupling for the case of symmetric quark matter with $\mu_{Q}=0$. At chemical
potentials below this critical value quark matter exists in a normal state,
while $\mu_{B}>\mu_{B}^{c}$ corresponds to 2SC phase. Larger diquark couplings
lead to stronger pairing among quarks and, consequently, yield smaller
$\mu_{B}^{c}$. At a certain value of the diquark coupling $\eta_{D}^{*}$ the
critical chemical potential vanishes, meaning that color superconductivity is
supported by the vacuum. Numerical analysis yields $\eta_{D}^{*}=0.78$. We
note that this result is not affected by the value of $\eta_{V}$ since vector
repulsion is irrelevant in the vacuum. In order to prevent formation of color
superconductivity in the vacuum we set the constraint $\eta_{D}<\eta_{D}^{*}$.
It is interesting to compare $\eta_{D}^{*}$ to the maximal value of the
diquark coupling for which color superconductivity in the vacuum is still
excluded. In the NJL model this is $({3}{m^{*}}/[2(m^{*}-m)]\simeq{3}/{2}$
Zablocki et al. (2010). Our value of $\eta_{D}^{*}$ is about twice smaller
than the NJL one. This difference is due to the adopted definition of
$\eta_{D}$. It can be exactly compensated by rescaling $\eta_{D}$ by the
factor ${G_{S0}}/{G_{PS0}}$. Indeed, in the vacuum $f_{{\bf k}fc}^{\pm}=0$ and
Eq. (II) at vanishing temperature, chemical potential and diquark pairing gap
yields

Figure 5: Baryonic chemical potential of the 2SC phase onset $\mu_{B}^{c}$ as
a function of the scaled diquark coupling $\eta_{D}$. The calculations are
performed for symmetric quark matter and $\eta_{V}=0$.

$\displaystyle\langle\overline{q}q\rangle=-6\sum_{f}\int\frac{d{\bf
k}}{(2\pi)^{3}}g_{\bf k}\frac{m^{*}}{\epsilon_{{\bf k}f}}.$ (42)

Here and all along the derivation of $\eta_{D}^{*}$ we omit the subscript
$``0"$ where it is applicable for the sake of shortening notations. Using this
relation we obtain by a straightforward calculation

$\displaystyle\frac{\partial^{2}\Omega}{\partial\Delta^{2}}\biggl{|}_{\Delta=0}$
$\displaystyle=$ $\displaystyle\frac{1}{2G_{D}}-4\sum_{f}\int\frac{d{\bf
k}}{(2\pi)^{3}}\frac{g_{\bf k}}{\epsilon_{{\bf k}f}}$ (43) $\displaystyle=$
$\displaystyle\frac{1}{2G_{D}}+\frac{2\langle{q}q\rangle}{3m^{*}}=0,$

This equation can be solved with respect to $G_{D}$. The corresponding
solution combined with $m^{*}=m-2G_{PS}\langle\overline{q}q\rangle$ and the
definition of $\eta_{D}$ yields

$\displaystyle\eta_{D}^{*}=\frac{3}{2}\frac{m^{*}}{m^{*}-m}\frac{G_{PS}}{G_{S}}.$
(44)

This expression coincides with the one of the NJL model only if
$G_{S}=G_{PS}$, which is violated by the vacuum solution of the present model.
At the same time, we would like to stress that the absolute value of
$G_{D}^{*}=\eta_{D}^{*}G_{S0}$ is about twice larger compared to the maximal
value of the diquark coupling in the NJL model. This provides the present
approach with a unique possibility to access the region of strong diquark
pairing that is interesting from the phenomenological point of view.

Figure 6: Effective quark mass $m^{*}$ (thick curves) and diquark pairing gap
$\Delta$ (thin curves) as functions of baryonic chemical potential $\mu_{B}$
calculated for electrically neutral $\beta$-equilibrated matter and
$\eta_{V}=0$.

##  VII Compact stars with quark cores

In order to model compact stars as astrophysical objects observable as
pulsars, we need the $T=0$ EoS under the conditions of electric neutrality and
$\beta$-equilibrium. The former is provided by a proper amount of non-
interacting electrons with the chemical potential $\mu_{e}$ and the negligible
mass $m_{e}=0.511$ MeV, while the latter is ensured by
$\mu_{e}=\mu_{d}-\mu_{u}=-\mu_{Q}$. The effective quark mass and the diquark
pairing gap calculated under these conditions are shown as functions of the
baryonic chemical potential in Fig. 6. At small $\mu_{B}$ chiral symmetry is
broken, which leads to high quark masses and the absence of color
superconductivity, as is indicated by a vanishing pairing gap. At a certain
value of the chemical potential $\mu_{B}$ the effective mass $m^{*}$ drops and
the diquark pairing gap $\Delta$ attains a finite value. The discontinuous
character of this change signals a first order phase transition. Further
increase of $\mu_{B}$ leads to a relatively slow decrease of effective mass
and a further growth of the pairing gap. At asymptotically high densities
holds $m^{*}\propto(\Delta\sum_{f}\mu_{f}^{*})^{-1}$ and
$\Delta\propto\sum_{f}\mu_{f}^{*3}$, see the Appendix B for the details. In
full agreement with the discussion of $\mu_{B}^{c}$, an increase of the
diquark coupling lowers the chemical potential of the chiral symmetry
restoration and color superconductivity onset. This picture holds for any
value of the vector coupling, which simply renormalizes $\mu_{B}$.

Having vector field, diquark pairing gap and chiral condensate self
consistently found from Eqs. (II) - (II), we can construct the EoS of
electrically neutral $\beta$-equilibrated quark matter. The hybrid EoS of
stellar matter is obtained by merging this quark matter EoS to the hadronic
one via the Maxwell construction corresponding to the first order phase
transition. The position of this transition is defined by requiring equal
pressures of quark and hadron phases at the same value of baryonic chemical
potential. In this case baryon number and energy densities experience a
discontinuous change along the transition. Its details are sensitive to the
values of vector and diquark couplings. In this work hadronic matter is
described by the DD2 EoS in its formulation and parametrization from Ref.
Typel et al. (2010). We use its version with hyperons DD2npY-T Shahrbaf et al.
(2022) that is supplemented with a crust EoS. This relativistic density
functional model introduces interaction between nucleons via meson exchange.
Besides reproducing the properties of normal nuclear matter, the DD2 EoS is in
a good agreement with the low density calculations of the chiral EFT approach
Krüger et al. (2013). The high density behavior of the DD2npY-T model can be
improved with the excluded volume correction Typel (2016). Within our
approach, however, this correction is neglected since at high densities
hadronic degrees of freedom are replaced by the quark ones. Fig. 7 shows the
hybrid quark-hadron EoS in the plane of energy density and pressure. Notably,
the onset of the 2SC phase in quark matter occurs at densities typical for the
hadron branch of the hybrid EoS. Therefore, the quark matter branch is already
color superconducting. Values of densities of the 2SC phase onset in quark
matter depend on the vector and diquark couplings and are given in the caption
of Fig. 7. As is seen from its upper panel, larger diquark couplings lead to
an earlier onset of quark matter. The width of the mixed phase region,
however, depends on $\eta_{D}$ only weakly. We would like to stress, that the
later the quark matter onset, the softer is its EoS. At vanishing $\eta_{V}$,
the quark part of the EoS is too soft and inconsistent with constraints
obtained within the multipolytrope analysis of the observational data of PSR
J1614+2230 Hebeler et al. (2013) and PSR J0740+6620 Miller et al. (2021). This
redundant softening can be compensated by the vector mean field repulsion in
quark matter. The corresponding effect is seen on the lower panel of Fig. 7.
It is remarkable, that the stiffening of the quark EoS can be so strong that
the corresponding pressure can sizably exceed the hadronic one. This yields a
positive feedback to the problem of reaching the two solar mass limit by
compact stars. The corresponding sets of the vector and diquark couplings lead
to an early onset of quark matter with the densities either lower or very
close to the one at which hyperons appear in the purely hadronic matter (white
circles on Fig. 7). As a result, within our approach the hyperon puzzle
Chatterjee and Vidaña (2016) is resolved by early onset of quark matter, which
already has been proposed as a possible solution, see Shahrbaf et al. (2020)
for a recent work and more literature cited therein. As a generic feature of
the present model we can conclude that the stiffness of the quark matter EoS
is controlled by the vector coupling, while the diquark one regulates the
quark matter onset density. Our hybrid EoS is in a good agreement with the
constraints from Ref. Hebeler et al. (2013). A similar constraint from Ref.
Miller et al. (2021) is fulfilled by the present approach only at high
densities. Some discrepancy at small densities is caused by the strong phase
transition from hadrons to quarks.

Figure 7: Hybrid EoS of cold electrically neutral $\beta$-equilibrated quark-
hadron matter in the plane of energy density $\varepsilon$ and pressure $p$.
Empty circles on the hadronic curves indicate the hyperon onset. The
calculations are performed for $\eta_{V}=0$ (upper panel) and $\eta_{V}\neq 0$
(lower panel). The nuclear matter constraints represented by the shaded areas
are discussed in the text. Arrows indicate the energy densities of the 2SC
phase onset in pure quark matter EoS, being a part of the hybrid EoS curve of
the corresponding color. These energy densities are 706 $\rm MeV~{}fm^{-3}$,
657 $\rm MeV~{}fm^{-3}$, 582 $\rm MeV~{}fm^{-3}$ and 546 $\rm MeV~{}fm^{-3}$
for the upper panel and 422 $\rm MeV~{}fm^{-3}$, 370 $\rm MeV~{}fm^{-3}$, 347
$\rm MeV~{}fm^{-3}$ and 313 $\rm MeV~{}fm^{-3}$ for the lower one.

This allows us to conclude about the limited ability of multipolytropic EoS to
catch a strong first order phase transition from hadronic to quark matter. It
is worth mentioning that a nonvanishing $\eta_{V}$ narrows the mixed phase
region when compared to the case of $\eta_{V}=0$.

The stiffness of an EoS can be quantified by the (squared) speed of sound
$c_{S}^{2}={dp}/{d\varepsilon}$. As is shown in the Appendix B, at high
densities and $\eta_{V}=0$, color superconducting quark matter is
characterized by $c_{S}^{2}\rightarrow{1}/{5}$, while a finite $\eta_{V}$
yields $c_{S}^{2}\rightarrow 1$. Fig. 8 shows this quantity as a function of
energy density. In the absence of vector repulsion and at densities typical
for the compact star interiors, the squared speed of sound of quark matter
varies within the range 0.3-0.35.

At vanishing vector coupling only the second of these scenarios can be
realized if quark matter onsets before two saturation densities provided by
$\eta_{D}>0.45$. However, this range of values of the diquark coupling is
inconsistent with the constraints on the stellar matter EoS from Refs. Hebeler
et al. (2013); Miller et al. (2021). This discrepancy can be resolved by
applying finite vector couplings. The lower panel of Fig. 8 indicates, that
the speed of sound increases when the vector coupling is increased, reflecting
a stiffening of the quark matter EoS. It is remarkable that for a given value
of $\eta_{V}$ the speed of sound of quark matter remains almost constant in
the density range typical for compact stars. For $\epsilon\leq 2000~{}{\rm
MeV~{}fm}^{-3}$ the corresponding variation does not exceed few percent. In
other words, our approach supports the constant speed of sound (CSS)
parameterization Alford et al. (2013); Zdunik and Haensel (2013), which is
widely applied to model quark matter, to classify compact stars with quark
cores Blaschke and Cierniak (2021); Miao et al. (2020) and to study the third
and fourth stellar families Paschalidis et al. (2018); Alford and Sedrakian
(2017); Li et al. (2020). It has been found recently that the EoS of color
superconducting nonlocal chiral quark models of the NJL type results in an
almost constant sound speed with $c_{S}^{2}$ in the narrow range between 0.45
and 0.54 Antić et al. (2021); Contrera et al. (2022).

We apply the family of developed EoS to model compact stars with quark cores
by solving the problem of relativistic hydrostatic equilibrium represented by
the TOV equations with the corresponding boundary conditions. Fig. 9 shows the
resulting relation between the stellar radius $\rm R$ and mass $\rm M$. Our
calculations are confronted with the constraint on the lower limit of the TOV
maximum mass given by the mass $2.01^{+0.04}_{+0.04}~{}\rm M_{\odot}$ measured
in a binary system of the pulsar PSR J0348+0432 and its white dwarf companion
Antoniadis et al. (2013). In the high and intermediate mass regions we also
utilize constraints from the Bayesian analysis of the observational data from
PSR J0740+6620 Riley et al. (2021); Miller et al. (2021) and PSR J0030+0451
Riley et al. (2019); Raaijmakers et al. (2019), respectively. A very important
information about the intermediate mass region of the $\rm M$ \- $\rm R$
diagram comes from the analysis of the gravitational wave signal produced by
the merger GW170817 Abbott et al. (2018). It constrains masses and radii of
the low and high mass objects in the binary as represented in Fig. 9 by the
orange and cyan shaded areas. The darker and lighter of them correspond to
$1\sigma$ and $2\sigma$ confidence intervals. Finally, based on the analysis
of the same gravitational wave signal we limit the star radius at $1.6~{}\rm
M_{\odot}$ from below by $R_{1.6}\geq 10.68$ km Bauswein et al. (2017) and at
$1.4~{}\rm M_{\odot}$ from above by $R_{1.4}\leq 13.6$ km Annala et al.
(2018).

Figure 8: Squared speed of sound $c_{S}^{2}$ of stellar matter as a function
of energy density $\varepsilon$ calculated with the quark-hadron EoS presented
in Fig. 7. Empty circles on the hadronic curves indicate the hyperon onset.

Regardless the values of $\eta_{D}$ and $\eta_{V}$, hybrid neutron stars with
quark cores that are obtained by a Maxwell construction are more compact
compared to their purely hadronic counterparts of the same mass. This effect
is especially attractive from the phenomenological point of view in the
typical stellar mass range of about $1.4~{}{\rm M}_{\odot}$, where the recent
observational data (in particular from GW170817 Abbott et al. (2018))
constrain the stellar radii to quite small values Capano et al. (2020);
Dietrich et al. (2020). Since our baseline hadronic EoS model DD2npY-T alone
would not fulfill these new radius constraints, this explains why the
phenomenologically interesting scenarios correspond to $\rm M_{onset}\lesssim
1.0~{}M_{\odot}$. The realization of such scenarios requires a relatively
large $\eta_{D}$ because that would lead to an early onset of deconfinement.
For the case of the three-flavor color superconducting NJL model see Klähn et
al. (2007, 2013). For the present model, this fact is illustrated in the upper
panel of Fig. 9, where an increase of the diquark coupling dramatically lowers
the onset mass of quark matter.

However, at vanishing vector coupling, the model does not support the
existence of stable hybrid stellar configurations with quark matter cores with
masses above $2\rm M_{\odot}$. Moreover, in this case stable hybrid star
branches the mass-radius diagram occur only if $\eta_{D}>0.43$.

We note that in this case, the remarkable feature of a third family of hybrid
stars is realised Gerlach (1968); Alvarez-Castillo et al. (2019); Blaschke et
al. (2020a) because the stable second family branch of pure neutron stars is
separated from the branch of stable hybrid stars by a sequence of unstable
stars that is characterised by $\partial\rm M/\partial\rm\varepsilon_{c}<0$
with $\varepsilon_{c}$ being energy density in the center of the star. The
existence of a third family is necessarily connected to the phenomenon of mass
twin stars Glendenning and Kettner (2000); Benic et al. (2014); Alvarez-
Castillo and Blaschke (2015); Zacchi et al. (2017); Alvarez-Castillo and
Blaschke (2017); Alford and Sedrakian (2017); Montana et al. (2019); Blaschke
et al. (2020b) the detection of which would be an observational proof of the
existence of a strong phase transition in dense neutron star matter.

The problem with a too low maximum mass can be overcome with adopting finite
values of the vector coupling. Increasing $\eta_{V}$ stabilizes quark branch
of the mass-radius diagram and causes a positive feedback to the problem of
reaching the two solar mass limit. However, the increase of $\eta_{V}$ moves
$\rm M_{onset}$ towards higher stellar masses, which requires simultaneous
increase of $\eta_{D}$ in order to keep $\rm M_{onset}<1.4~{}M_{\odot}$. We
note that this requirement is set in order to obtain the necessary softening
of the hybrid EoS for fulfilling the tidal deformability constraint formulated
for $1.4~{}M_{\odot}$ stars, which is not fulfilled by the DD2npY - T EoS
itself. For example, at $\eta_{D}=0.530$ the maximum mass $\rm M_{max}$ of the
stable hybrid star sequence reaches $2.01~{}{\rm M_{\odot}}$ for
$\eta_{V}=0.172$, while marginal agreement with the PSR J0740+6620
observational data is reached at $\eta_{V}=0.177$. Further increase of the
vector coupling allows the present model to fulfil the above-mentioned
astrophysical constraints. It is also interesting to note that at a certain
value of $\eta_{V}$ defined by the value of $\eta_{D}$, the hybrid star branch
of the mass-radius curve intersects the stable part of the branch
corresponding to purely hadronic EoS models of neutron star matter. For
$\eta_{D}=0.530$ this happens at $\eta_{V}=0.216$. In other words, at
sufficiently large vector and diquark couplings, the existence of quark matter
cores in the compact star interiors allows them to break through to the region
of the mass-radius diagram inaccessible for purely hadronic stars due to the
softening of their EoS caused by the appearance of hyperons. This possibility
is important for compact star phenomenology since it might provide a hybrid
neutron star explanation of the nature of the $2.6~{}\rm M_{\odot}$ companion
of the $23~{}\rm M_{\odot}$ black hole recently reported in the merger
GW190814 Abbott et al. (2020). We discuss this issue below.

High values of the vector and diquark couplings also allow the present model
to achieve simultaneous agreement with the high and intermediate mass
constraints coming from the NICER data on PSR J0740+6620 and from the
gravitational wave signal of GW170817, which is a challenging task for most of
the other approaches.

Figure 9: Mass-radius relation of hybrid neutron stars with the quark-hadron
EoS presented on Fig. 7. An empty circle on the hadronic curves indicates the
hyperon onset. The blue filled circles represent the special points with the
mass $\rm M_{SP}$ found according to the fitting procedure described in the
text. The astrophysical constraints depicted by the colored bends and shaded
areas are discussed in the text.

Fig. 9 demonstrates a special aspect of star sequences in the mass-radius
diagram that occurs only for hybrid stars, the so-called ”special point”.
Strictly speaking, it is not a point but a narrow region of intersection of
hybrid star sequences, which was first reported and analysed in Ref. Yudin et
al. (2014). A thorough of its invariance properties has been given recently in
Refs. Cierniak and Blaschke (2020); Blaschke and Cierniak (2021); Cierniak and
Blaschke (2021, 2022). The remarkable feature of the insensitivity of the
special point location to the details of hadronic EoS and the construction of
the phase transition, makes it a universal tool to study the properties of
quark cores in compact stars. Within the CSS parameterization this special
point is insensitive to constant bag pressure controlling the quark matter
onset density. In our approach the same phenomenological role is played by the
diquark coupling. This explains the occurrence of the special point when
varying just $\eta_{D}$ in the upper panel of Fig. 9. Its mass $M_{\rm SP}$
was found to obey empirical relation

$\displaystyle\rm M_{max}=M_{SP}+\delta(M_{SP}-M_{onset})$ (45)

with $\delta$ being a constant Blaschke and Cierniak (2021). The value of this
parameter depend on a particular choice of the quark matter EoS. For the CSS
parameterization holds $\delta=0.1$ Blaschke and Cierniak (2021). We
independently find $\rm M_{SP}$ and $\delta$ by fitting Eq. (45) to $\rm
M_{max}$ and $\rm M_{onset}$ extracted from the mass-radius diagrams shown on
Fig. 9. Note that for the $\eta_{V}=0$ case we consider only the mass-radius
curves supporting stable quark branches. In this case $\rm
M_{SP}=1.31M_{\odot}$ and $\delta=0.06$. For the case of finite vector
coupling shown on the lower panel of Fig. 9, we obtain $\rm
M_{SP}=1.97M_{\odot}$ and $\delta=0.24$. As is seen from Fig. 9, $\rm M_{SP}$
perfectly fits the special point mass at vanishing and finite $\eta_{V}$.

The gravitational wave signal frrom the inspiral phase of the binary neutron
star merger GW170817 is another source of information about the EoS for
neutron star matter. The most important quantity extracted form the
observational data is the dimensionless tidal deformability
$\Lambda=\frac{2}{3}k_{2}C^{2}$ expressed through the second Love number
$k_{2}$ and the stellar compactness $C=\frac{\rm M}{\rm R}$. For a compact
star with the mass $\rm M=1.4~{}M_{\odot}$ it is constrained to
$\Lambda_{1.4}=190^{+390}_{-120}$ Abbott et al. (2018). At vanishing vector
coupling, $\eta_{V}=0$, the EoS of quark matter does not support this value of
the stellar mass. Therefore, we analyse $\Lambda$ only at finite vector
coupling. The dependence of the dimensionless tidal deformability on the
stellar mass is shown on the top panel of Fig. 10. The observational data
constraint on $\Lambda$ is fulfilled only for sufficiently high diquark
couplings providing a low onset mass of quark matter and thus a sufficient
compactness of the stellar configurations. This constraint is fulfilled by the
curves (0.235, 0.540) and (0.265, 0.555) yielding $\Lambda_{1.4}=395$ and
$\Lambda_{1.4}=325$, respectively. The best agreement is achieved for the same
model, which is the most accurate with respect to the discussed observational
constraint on the mass-radius diagram.

The masses of two components $\rm M_{1}$ and $\rm M_{2}$ of the event GW170817
are connected by the chirp mass Peters and Mathews (1963)

$\displaystyle\mathcal{M}=\frac{(\rm M_{1}M_{2})^{\frac{3}{5}}}{(\rm
M_{1}+M_{2})^{\frac{1}{5}}},$ (46)

which is an important characteristics of the merger. By convention $\rm
M_{1}\geq M_{2}$. The best measured combination of the masses from GW170817
yields $\mathcal{M}=1.188^{+0.004}_{-0.002}~{}M_{\odot}$ Abbott et al. (2017),
which is used below. We also rely on the results of the recent analysis
imposing a common EoS of two components Abbott et al. (2018). Compared to the
case of independently imposed EoS Abbott et al. (2019), it gives a more
stringent constraint on the tidal deformabilities of two components
$\Lambda_{1}$ and $\Lambda_{2}$. For each set of $\rm M_{1}$ and $\rm M_{2}$
providing the chosen value of $\mathcal{M}$, the tidal deformabities of two
components were calculated through EoS-insensitive relations Yagi and Yunes
(2013). The results are demonstrated in the lower panel of Fig. 10. Increasing
the diquark coupling lowers the quark matter onset mass, leading in turn to
more compact stellar configurations and smaller tidal deformabilities.
However, an increase of the vector coupling stiffens the quark matter EoS and
causes the tidal deformability of the star to grow. Two curves fulfilling the
constraint on $\Lambda_{1.4}$, namely the ones for the parameter sets (0.235,
0.540) and (0.265, 0.555), are well inside the 90% confidence level and lie
not so far from the 50% one.

Figure 10: Dimensionless tidal deformability $\Lambda$ as a function of
stellar mass $\rm M$ (upper panel) and the same quantity of the low mass
component of the NS-NS merger $\Lambda_{2}$ as a function of the corresponding
parameter $\Lambda_{1}$ of the high mass one (lower panel). The calculations
are performed for the quark-hadron EoSs presented on the lower panel of Fig.
8. The $\Lambda_{1}-\Lambda_{2}$ curve $(0.180,0.505)$ is not shown since
within the shown region it coincides with the purely hadronic one. Colored
filled circles on the lower panel represent the configurations with equal
masses of two components $\rm M_{1}=M_{2}=1.3646~{}M_{\odot}$. Dark and light
green shaded areas represent the regions falling into the $50~{}\%$ and
$90~{}\%$ confidence levels.

The 2.5 - 2.67 $\rm M_{\odot}$ compact object from the merger GW190814 is
either the lightest black hole or the heaviest neutron star ever observed
Abbott et al. (2020). Despite the fact that both of these possibilities are
neither confirmed nor falsified, the latter one raises the question about
stability of ultra heavy neutron stars. As was mentioned before, sufficiently
high values of the vector and diquark couplings stabilize quark matter in the
interiors of heavy compact stars. However, the values of $\eta_{V}$ and
$\eta_{D}$ needed to provide stability of the 2.6 $M_{\odot}$ stars within our
approach prevent the Maxwell construction of the quark-hadron transition due
to absence of crossing of the quark matter EoS with the DD2npY - T hadronic
one in the $\mu_{B}-p$ plane. It is worth mentioning, that this issue itself
might be related to the phase transition construction, i.e. the Glendenning
Glendenning (1992), pasta Heiselberg et al. (1993); Yasutake et al. (2014);
Maslov et al. (2019) or interpolation Masuda et al. (2013); Alvarez-Castillo
et al. (2014); Baym et al. (2018); Blaschke and Alvarez-Castillo (2020);
Abgaryan et al. (2018); Ayriyan et al. (2021); Ivanytskyi and Blaschke (2022)
constructions. The corresponding analysis lays out of the scope of the present
work. Therefore, we use as a stiffer hadronic EoS the DD2p40 one being a
version of the DD2 EoS with the excluded volume correction Typel (2016) that
allows the Maxwell construction at high values of $\eta_{V}$ and $\eta_{D}$.
The mass-radius relation obtained with the corresponding hybrid EoS is shown
on Fig. 11. Remarkably, a proper choice of the vector and diquark couplings
allows our model not only to reach 2.6 $M_{\odot}$ stellar mass, but also
makes it consistent with the stellar radius and tidal deformability
constraints from the GW170817 merger Abbott et al. (2018). For example, at
$\eta_{V}=0.425$ and $\eta_{D}=0.575$ our approach yields $\rm
M_{max}=2.64~{}M_{\odot}$, $\rm R_{1.4}=12.6$ km and $\Lambda_{1.4}=543$. It
is also worth mentioning, that these values of the vector and diquark
couplings lead to a rather early onset of quark matter. The corresponding
values of the baryon density are $n_{B}=0.17~{}fm^{-3}$ on the hadron boundary
of the mixed phase and $n_{B}=0.27~{}fm^{-3}$ on the quark one. At the
considered values of the coupling constants and for densities typical for
stellar interiors, the squared speed of sound of quark matter varies within
the interval $c_{S}^{2}=0.63-0.8$.

As is seen from Fig. 11, the set of hybrid EoS selected has a special point
with the mass $\rm M_{SP}=2.56~{}M_{\odot}$ perfectly fitting Eq. (45) with
$\delta=0.035$. This value of $\delta$ differs from those found for the
special points shown on Fig. 9. This supports our conclusion about the model
dependence of the scaling (45), which itself is universal.

Figure 11:  The same as in Fig. 9 but with the DD2p40 model for hadron part of
hybrid EoS and with different values of $\eta_{V}$ and $\eta_{D}$ given in the
legend. The upper colored band represents the range of stellar masses reported
in Ref. Abbott et al. (2020).

##  VIII Conclusions

We have studied a relativistic density functional approach to quark matter,
which i) mimics quark confinement by a rapid growth of the quark self-energy
in the confining region, ii) respects chiral symmetry of strong interaction
Lagrangian model and iii) can be interpreted as a chiral quark model with
density dependent coupling constants. In addition to the vector repulsion
channel, we have introduced the diquark pairing that was not studied before
within the relativistic density functional approaches to quark matter.

In order to model NS with quark cores we apply the present density functional
approach to Maxwell-construct a hybrid EoS made of a hadronic one and the new
color superconducting quark matter EoS with the vector repulsion. The results
of this modeling are confronted with the recent constraints on the NS mass-
radius relation and tidal deformability from multi-messenger astronomy. As a
general feature, we want to stress that color superconductivity lowers the
onset density of quark matter and provides the present approach with the
possibility to describe hybrid star sequences with a maximum mass above 2.6
$\rm M_{\odot}$, while being sufficiently compact at 1.4 $\rm M_{\odot}$,
consistent with the mentioned constraints. The approach has still a potential
for further improvement, e.g., by adjusting the values of the vector and
diquark couplings and by modifying the hadron-to-quark matter transition
construction.

## Acknowledgements

This work has been supported by the Polish National Science Centre (NCN) under
grant No. 2019/33/B/ST9/03059. It is part of a project that has received
funding from the European Union’s Horizon 2020 research and innovation program
under grant agreement STRONG – 2020 - No 824093.

## APPENDIX A

Here we consider the confinement mechanism from Ref. Fowler et al. (1981) in
order to demonstrate how it relates to the ones from the present work and from
Ref. Kaltenborn et al. (2017). For this we start with the zero temperature
pressure of the gas of quark quasi particles with medium dependent mass
$m^{*}$ and chemical potentials $\mu_{f}$

$\displaystyle p=2\sum_{f,c}\int\frac{d{\bf
k}}{(2\pi)^{3}}\left[\epsilon_{{\bf k}f}g_{\bf k}+(\mu_{f}-\epsilon_{{\bf
k}f})f_{{\bf k}f}\right]-\mathcal{U}+n_{s}\frac{\partial\mathcal{U}}{\partial
n_{s}},$

where $f_{{\bf k}f}=\theta(\mu_{f}-\epsilon_{{\bf k}f})$, $\mathcal{U}$ is the
quasiparticle interaction potential, which is a function of their scalar
density $n_{s}$ discussed below. The remaining notations coincide with the
ones adopted in the main text of this work. We note that the last two terms in
Eq. (APPENDIX A) are equivalent to
$-\mathcal{U}_{MF}+\langle\overline{q}q\rangle\Sigma_{MF}$ in Eq. (8). At the
moment we do not discuss the form of $\mathcal{U}$ in order to keep the
analysis general. Contrary to the present consideration, Ref. Fowler et al.
(1981) neglects the interaction potential and the zero point energy term that
represents the contribution of fermionic zero mode.

The condition of the pressure stationarity $\partial p/\partial m^{*}=0$ is
equivalent to Eq. (II). The definition of the scalar density

$\displaystyle n_{s}=2\sum_{f,c}\int\frac{d{\bf k}}{(2\pi)^{3}}\left[f_{{\bf
k}f}-g_{\bf k}\right]\frac{m^{*}}{\epsilon_{{\bf k}f}}.$ (A.2)

allows us to write this condition in the form $\partial m^{*}/\partial
n_{s}=\partial^{2}\mathcal{U}/\partial n_{s}^{2}$. It can be integrated in
order to obtain the mass gap equation

$\displaystyle m^{*}=m+\frac{\partial\mathcal{U}}{\partial n_{s}},$ (A.3)

where the current quark mass $m$ appears as an integration constant. Eq. (A.3)
along with the definition (A.2) should be solved selfconsistently in order to
find $m^{*}$ for a given $\mathcal{U}$. Once this is done, the quark number
density $n$ and energy density $\varepsilon$ can be found using the
thermodynamic identities discussed in Section II

$\displaystyle n$ $\displaystyle=$ $\displaystyle 2\sum_{f,c}\int\frac{d{\bf
k}}{(2\pi)^{3}}f_{{\bf k}f},$ (A.4) $\displaystyle\varepsilon$
$\displaystyle=$ $\displaystyle 2\sum_{f,c}\int\frac{d{\bf
k}}{(2\pi)^{3}}\left[f_{{\bf k}f}-g_{\bf k}\right]\epsilon_{{\bf
k}f}+\mathcal{U}-n_{s}\frac{\partial\mathcal{U}}{\partial n_{s}}.$ (A.5)

Eq. (A.5) coincides with Eq. (1) of Ref. Fowler et al. (1981) if the
contributions of fermionic zero mode and interaction potential are neglected.
Let us consider the regime of vanishing quark number density. In this case
$m^{*}$ is large and the low-momentum states with $\epsilon_{{\bf k}f}\simeq
m^{*}$ are the most important. This allows us to approximate

$\displaystyle n_{s}$ $\displaystyle=$ $\displaystyle
n-2\sum_{f,c}\int\frac{d{\bf k}}{(2\pi)^{3}}g_{\bf k},$ (A.6)
$\displaystyle\varepsilon$ $\displaystyle=$ $\displaystyle
m^{*}n_{s}+\mathcal{U}-n_{s}\frac{\partial\mathcal{U}}{\partial n_{s}}.$ (A.7)

Adopting the assumption of Ref. Fowler et al. (1981) that $\varepsilon$ at
vanishing density converges to some constant value $B$, we conclude from the
last expression that

$\displaystyle
m^{*}=\frac{B-\mathcal{U}}{n_{s}}+\frac{\partial\mathcal{U}}{\partial n_{s}}.$
(A.8)

This result demonstrates that at vanishing $\mathcal{U}$ the low density
scaling $m^{*}=B/n$ from Ref. Fowler et al. (1981) is respected only at
$n_{s}=n$ provided by the no-sea approximation that neglects the contribution
of fermionic zero mode ($g_{\bf k}=0$). On the other hand, setting
$\mathcal{U}=0$ at any $m^{*}\neq m$ does not fulfill the mass gap equation
(A.3). If it is respected, then Eq. (A.7) turns to
$\varepsilon=mn_{s}+\mathcal{U}$, which along with $\varepsilon=B$ makes Eqs.
(A.3) and (A.8) coincide. For the SFM potential
$\mathcal{U}_{SFM}=D_{0}n_{s}^{2/3}$ this yields
$m^{*}_{SFM}=m+(2/3)D_{0}n_{s}^{-1/3}$, which at small densities and $g_{\bf
k}=0$ gives $m^{*}_{SFM}\propto n^{-1/3}$ Horowitz et al. (1985); Röpke et al.
(1986); Kaltenborn et al. (2017). Thus, the low density scaling $m^{*}\propto
n^{-1}$ requires the no-sea approximation and violates the mass gap equation
that would be required for the pressure stationarity.

## APPENDIX B

Here we analyze the high density asymptotic of the 2SC phase at vanishing
temperature. For definiteness we consider $\mu_{f}^{*}>0$. In this regime
$\mu_{f}^{*}\gg\Lambda\gg m_{f}^{*}$ and quark masses can be neglected as well
as zero point terms. In addition, only quarks ($a=+$) with $f_{{\bf
k}fc}^{+}=\theta(\mu_{f}^{*}-|{\bf k}|)$ contribute to the thermodynamic
quantities. Then, the pairing gap equation becomes

$\displaystyle\Delta\simeq-4G_{D}\sum_{f,c}\int\frac{d{\bf
k}}{(2\pi)^{3}}f_{{\bf k}fc}^{+}\frac{\Delta_{c}}{\epsilon_{{\bf
k}fc}^{+}}.\,$ (B.1)

At $\Delta\sim\mu_{f}^{*}$ or at $\Delta\ll\mu_{f}^{*}$, the right hand side
of this equation diverges as $\sim\Delta\sum_{f}\mu_{f}^{*2}$, which is much
faster than $\Delta$ corresponding to the left hand side. Therefore,
$\Delta\gg\mu_{f}^{*}$. In this case the upper edge of the filled single
particle states lies below the Fermi surface, $\epsilon_{{\bf
k}fc}^{+}\simeq-\Delta$, and Eq. (B.1) becomes

$\displaystyle\Delta\simeq 4G_{D}\sum_{f,c}\int\frac{d{\bf
k}}{(2\pi)^{3}}f_{{\bf
k}fc}^{+}\frac{\Delta_{c}}{\Delta}=\frac{4G_{D}}{3\pi^{2}}\sum_{f}\mu_{f}^{*3}.$
(B.2)

Similarly, the equation for the vector field simplifies to

$\displaystyle\omega\simeq-4G_{V}\sum_{f,c}\int\frac{d{\bf
k}}{(2\pi)^{3}}f_{{\bf k}fc}^{+}\left(2\Delta_{c}\delta(\epsilon_{{\bf
k}fb}^{+})+\frac{\epsilon_{{\bf k}fb}^{+}}{\epsilon_{{\bf k}fc}^{+}}\right).$
(B.3)

For paired quarks, the bracket in this expression can be approximated by
$2\Delta\delta(\epsilon_{{\bf k}fb}^{+})$ since $\frac{\epsilon_{{\bf
k}fb}^{+}}{\epsilon_{{\bf k}fc}^{+}}\ll 1$, while $f_{{\bf k}fc}^{+}$ can be
replaced by one half due to the Dirac delta function. The corresponding
momentum integral behaves as $\sim\Delta\sum_{f}\mu_{f}^{*2}$. The
contribution of unpaired blue quarks can be neglected since it behaves as
$\sim\sum_{f}\mu_{f}^{*3}$. Thus, with $\epsilon_{{\bf k}fb}^{+}\simeq|{\bf
k}|-\mu_{f}^{*}$ the vector field becomes

$\displaystyle\omega\simeq-4G_{V}\Delta\sum_{f,c=r,g}\int\frac{d{\bf
k}}{(2\pi)^{3}}\delta(\epsilon_{{\bf
k}fb}^{+})=-\frac{4G_{V}\Delta}{\pi^{2}}\sum_{f}\mu_{f}^{*2}.$ (B.4)

The quark number density is

$\displaystyle\langle
q^{+}q\rangle=-\frac{\omega}{2G_{V}}\simeq\frac{2\Delta}{\pi^{2}}\sum_{f}\mu_{f}^{*2}.$
(B.5)

This gives a direct access to number density of a given quark flavor $\langle
f^{+}f\rangle\simeq 2\Delta\mu_{f}^{*2}/\pi^{2}$.

In order to find the asymptotic value of the effective quark mass we also
consider the chiral condensate. Applying the simplifications adopted before
and $\epsilon_{{\bf k}f}\simeq|{\bf k}|$ one gets

$\displaystyle\langle\overline{q}q\rangle\simeq 4\Delta
m^{*}\sum_{f}\int\frac{d{\bf k}}{(2\pi)^{3}}\frac{\delta(\epsilon_{{\bf
k}fb}^{+})}{|{\bf k}|}=\frac{2\Delta m^{*}}{\pi^{2}}\sum_{f}\mu_{f}^{*}.$
(B.6)

Using definition of the effective quark mass, we obtain

$\displaystyle
m^{*}=m-2G_{PS}\langle\overline{q}q\rangle\simeq\frac{m}{1+\frac{4G_{PS}\Delta}{\pi^{2}}\sum\limits_{f}\mu_{f}^{*}}.$
(B.7)

The quark pressure can be found using Eq. (19) as
$p_{q}=-\Omega_{q}+\Omega_{q0}$, with $\Omega_{q0}$ being quark thermodynamic
potential in the vacuum. For this we integrate by parts, neglect
$\Omega_{q0}$, zero point and antiquark terms and arrive to

$\displaystyle p_{q}\simeq 2\sum_{f,c}\int\frac{d{\bf
k}}{(2\pi)^{3}}f^{+}_{{\bf k}fc}\left(2\Delta_{c}\delta(\epsilon_{{\bf
k}fb}^{+})+\frac{\epsilon_{{\bf k}fb}^{+}}{\epsilon_{{\bf
k}fc}^{+}}\right)\frac{{\bf k}^{2}}{3\epsilon_{{\bf k}f}}.$ (B.8)

The contribution of paired quarks can be found similar to the case of
$\omega$. It is
$\frac{2\Delta}{3\pi^{2}}\sum_{f}\mu_{f}^{*3}\simeq\frac{\Delta^{2}}{2G_{D}}$,
where Eq. (B.2) was used. The blue unpaired quarks contribute to $p_{q}$ as
$\sim\frac{1}{12\pi^{2}}\sum_{f}\mu_{f}^{*4}$, which can be neglected. Thus,
the total pressure becomes

$\displaystyle p\simeq
p_{q}+\frac{\omega^{2}}{4G_{V}}-\frac{\Delta^{2}}{4G_{D}}\simeq\frac{\omega^{2}}{4G_{V}}+\frac{\Delta^{2}}{4G_{D}}.$
(B.9)

The energy density is expressed using the corresponding thermodynamic identity
and the definition of the effective chemical potential

$\displaystyle\varepsilon=\sum_{f}\langle
f^{+}f\rangle\mu_{f}-p\simeq\frac{\omega^{2}}{4G_{V}}+\frac{5\Delta^{2}}{4G_{D}}.$
(B.10)

At finite vector coupling $\omega\gg\Delta$ and $\varepsilon\simeq
p\simeq\frac{\omega^{2}}{4G_{V}}$, while at $G_{V}=0$, one gets $\omega=0$ and
$\varepsilon\simeq 5p\simeq\frac{5\Delta^{2}}{4G_{D}}$. This yields the
asymptotic value of the squared speed of sound

$\displaystyle
c_{S}^{2}=\frac{dp}{d\varepsilon}\simeq\left\\{\begin{array}[]{l}1,\quad
G_{V}\neq 0\\\ \frac{1}{5},\quad G_{V}=0\end{array}\right..$ (B.13)

Now we focus on the electrically neutral case of $\beta$-equilibrium at
$G_{V}\neq 0$. Using above mentioned expressions for number density of $u$ and
$d$ quarks along with the electron number density
$n_{e}=\mu_{e}^{3}/(3\pi^{2})$ the corresponding condition becomes

$\displaystyle\frac{2}{3}\frac{2\Delta\mu_{u}^{*2}}{\pi^{2}}-\frac{1}{3}\frac{2\Delta(\mu_{u}^{*}+\mu_{e})^{2}}{\pi^{2}}-\frac{\mu_{e}^{3}}{3\pi^{2}}=0.$
(B.14)

where $\mu_{d}^{*}=\mu_{u}^{*}+\mu_{e}$ is used in order to express effective
chemical potential of $d$ quarks. The electron term in Eq. (B.14) is
negligibly small compared to the quark ones since
$\mu_{e}\sim\mu_{f}^{*}\ll\Delta$. With this we find
$\mu_{e}\simeq\left(2^{1/2}-1\right)\mu_{u}^{*}$ and $\mu_{d}^{*}\simeq
2^{1/2}\mu_{u}^{*}$. Inserting this $\mu_{d}^{*}$ to Eqs. (B.2) and (B.4) we
find vector field $\omega\sim\mu_{u}^{*5}$ and effective chemical potential of
$u$ quarks

$\displaystyle\mu_{u}^{*}\simeq\mu_{u}-\frac{16G_{D}G_{V}}{3\pi^{4}}\mu_{u}^{*5}.$
(B.15)

Two terms on the right hand side of this equation diverge faster then
$\mu_{u}^{*}$ on the left hand side, which can be neglected. This yields
$\mu_{u}^{*}\sim\mu_{u}^{1/5}$. We also approximate chemical potential of $u$
quarks as $\mu_{u}=\mu_{B}/3-2\mu_{e}/3\simeq\mu_{B}/3$, which is valid since
$\mu_{e}\sim{\mu_{u}}^{1/5}\ll\mu_{u}$. Thus

$\displaystyle\mu_{u}^{*}\simeq\left(\frac{\pi^{4}\mu_{B}}{16G_{D}G_{V}}\right)^{\frac{1}{5}},$
(B.16)

With this equation and $\mu_{d}^{*}=2^{1/2}\mu_{u}^{*}$ asymptotic expressions
for pairing gap and quark number density become

$\displaystyle\Delta$ $\displaystyle\simeq$
$\displaystyle\frac{4G_{D}}{3\pi^{2}}\left(1+2^{\frac{3}{2}}\right)\left(\frac{\pi^{4}\mu_{B}}{16G_{D}G_{V}}\right)^{\frac{3}{5}},$
(B.17) $\displaystyle\langle q^{+}q\rangle$ $\displaystyle\simeq$
$\displaystyle\frac{1+2^{\frac{3}{2}}}{2G_{V}}\mu_{B}.$ (B.18)

## References

  * Bazavov et al. (2019) A. Bazavov et al. (HotQCD), Phys. Lett. B 795, 15 (2019), eprint 1812.08235.
  * Glendenning (1992) N. K. Glendenning, Phys. Rev. D 46, 1274 (1992).
  * Heiselberg et al. (1993) H. Heiselberg, C. J. Pethick, and E. F. Staubo, Phys. Rev. Lett. 70, 1355 (1993).
  * Yasutake et al. (2014) N. Yasutake, R. Lastowiecki, S. Benic, D. Blaschke, T. Maruyama, and T. Tatsumi, Phys. Rev. C 89, 065803 (2014), eprint 1403.7492.
  * Maslov et al. (2019) K. Maslov, N. Yasutake, A. Ayriyan, D. Blaschke, H. Grigorian, T. Maruyama, T. Tatsumi, and D. N. Voskresensky, Phys. Rev. C 100, 025802 (2019), eprint 1812.11889.
  * Masuda et al. (2013) K. Masuda, T. Hatsuda, and T. Takatsuka, PTEP 2013, 073D01 (2013), eprint 1212.6803.
  * Alvarez-Castillo et al. (2014) D. E. Alvarez-Castillo, S. Benić, D. Blaschke, and R. Łastowiecki, Acta Phys. Polon. Supp. 7, 203 (2014), eprint 1311.5112.
  * Baym et al. (2018) G. Baym, T. Hatsuda, T. Kojo, P. D. Powell, Y. Song, and T. Takatsuka, Rept. Prog. Phys. 81, 056902 (2018), eprint 1707.04966.
  * Blaschke and Alvarez-Castillo (2020) D. Blaschke and D. Alvarez-Castillo, Eur. Phys. J. A 56, 124 (2020), eprint 1807.03258.
  * Abgaryan et al. (2018) V. Abgaryan, D. Alvarez-Castillo, A. Ayriyan, D. Blaschke, and H. Grigorian, Universe 4, 94 (2018), eprint 1807.08034.
  * Ayriyan et al. (2021) A. Ayriyan, D. Blaschke, A. G. Grunfeld, D. Alvarez-Castillo, H. Grigorian, and V. Abgaryan, Eur. Phys. J. A 57, 318 (2021), eprint 2102.13485.
  * Ivanytskyi and Blaschke (2022) O. Ivanytskyi and D. Blaschke (2022), eprint 2205.03455.
  * Kaltenborn et al. (2017) M. A. R. Kaltenborn, N.-U. F. Bastian, and D. B. Blaschke, Phys. Rev. D 96, 056024 (2017), eprint 1701.04400.
  * Horowitz et al. (1985) C. J. Horowitz, E. J. Moniz, and J. W. Negele, Phys. Rev. D 31, 1689 (1985).
  * Röpke et al. (1986) G. Röpke, D. Blaschke, and H. Schulz, Phys. Rev. D 34, 3499 (1986).
  * Blaschke et al. (2022) D. Blaschke, O. Ivanytskyi, and M. Shahrbaf, _Quark deconfinement in compact stars through sexaquark condensation_ (2022), eprint 2202.05061.
  * Klähn et al. (2007) T. Klähn, D. Blaschke, F. Sandin, C. Fuchs, A. Faessler, H. Grigorian, G. Röpke, and J. Trümper, Phys. Lett. B 654, 170 (2007), eprint nucl-th/0609067.
  * Klähn et al. (2013) T. Klähn, R. Łastowiecki, and D. B. Blaschke, Phys. Rev. D88, 085001 (2013), eprint 1307.6996.
  * Blaschke et al. (2014) D. Blaschke, M. Buballa, A. Dubinin, G. Röpke, and D. Zablocki, Annals Phys. 348, 228 (2014), eprint 1305.3907.
  * Cooper (1956) L. N. Cooper, Phys. Rev. 104, 1189 (1956).
  * Klevansky (1992) S. P. Klevansky, Rev. Mod. Phys. 64, 649 (1992).
  * Hüfner et al. (1994) J. Hüfner, S. P. Klevansky, P. Zhuang, and H. Voss, Annals Phys. 234, 225 (1994).
  * Zhuang et al. (1994) P. Zhuang, J. Hüfner, and S. P. Klevansky, Nucl. Phys. A 576, 525 (1994).
  * Buballa (2005) M. Buballa, Phys. Rept. 407, 205 (2005), eprint hep-ph/0402234.
  * Gell-Mann et al. (1968) M. Gell-Mann, R. J. Oakes, and B. Renner, Phys. Rev. 175, 2195 (1968).
  * Fowler et al. (1981) G. N. Fowler, S. Raha, and R. M. Weiner, Z. Phys. C 9, 271 (1981).
  * Gross and Wilczek (1973) D. J. Gross and F. Wilczek, Phys. Rev. Lett. 30, 1343 (1973).
  * Politzer (1973) H. D. Politzer, Phys. Rev. Lett. 30, 1346 (1973).
  * Deur et al. (2016) A. Deur, S. J. Brodsky, and G. F. de Teramond, Prog. Part. Nucl. Phys. 90, 1 (2016), eprint 1604.08082.
  * Ratti et al. (2006) C. Ratti, M. A. Thaler, and W. Weise, Phys. Rev. D 73, 014019 (2006), eprint hep-ph/0506234.
  * Grigorian (2007) H. Grigorian, Phys. Part. Nucl. Lett. 4, 223 (2007), eprint hep-ph/0602238.
  * Jamin (2002) M. Jamin, Phys. Lett. B 538, 71 (2002), eprint hep-ph/0201174.
  * Tanabashi et al. (2018) M. Tanabashi, K. Hagiwara, K. Hikasa, K. Nakamura, Y. Sumino, F. Takahashi, J. Tanaka, K. Agashe, G. Aielli, C. Amsler, et al. (Particle Data Group), Phys. Rev. D 98, 030001 (2018), URL https://link.aps.org/doi/10.1103/PhysRevD.98.030001.
  * Plumari et al. (2011) S. Plumari, W. M. Alberico, V. Greco, and C. Ratti, Phys. Rev. D 84, 094004 (2011), eprint 1103.5611.
  * Bozek et al. (1998) P. Bozek, Y. B. He, and J. Hufner, Phys. Rev. C 57, 3263 (1998), eprint nucl-th/9802021.
  * Zablocki et al. (2010) D. S. Zablocki, D. B. Blaschke, R. Anglani, and Yu. L. Kalinovsky, Acta Phys. Polon. Supp. 3, 771 (2010), eprint 0912.4929.
  * Typel et al. (2010) S. Typel, G. Röpke, T. Klähn, D. Blaschke, and H. H. Wolter, Phys. Rev. C 81, 015803 (2010), eprint 0908.2344.
  * Shahrbaf et al. (2022) M. Shahrbaf, D. Blaschke, S. Typel, G. R. Farrar, and D. E. Alvarez-Castillo, Phys. Rev. D 105, 103005 (2022), eprint 2202.00652.
  * Krüger et al. (2013) T. Krüger, I. Tews, K. Hebeler, and A. Schwenk, Phys. Rev. C 88, 025802 (2013), eprint 1304.2212.
  * Typel (2016) S. Typel, Eur. Phys. J. A 52, 16 (2016).
  * Hebeler et al. (2013) K. Hebeler, J. M. Lattimer, C. J. Pethick, and A. Schwenk, Astrophys. J. 773, 11 (2013), eprint 1303.4662.
  * Miller et al. (2021) M. C. Miller et al., Astrophys. J. Lett. 918, L28 (2021), eprint 2105.06979.
  * Chatterjee and Vidaña (2016) D. Chatterjee and I. Vidaña, Eur. Phys. J. A 52, 29 (2016), eprint 1510.06306.
  * Shahrbaf et al. (2020) M. Shahrbaf, D. Blaschke, A. G. Grunfeld, and H. R. Moshfegh, Phys. Rev. C 101, 025807 (2020), eprint 1908.04740.
  * Alford et al. (2013) M. G. Alford, S. Han, and M. Prakash, Phys. Rev. D 88, 083013 (2013), eprint 1302.4732.
  * Zdunik and Haensel (2013) J. L. Zdunik and P. Haensel, Astron. Astrophys. 551, A61 (2013), eprint 1211.1231.
  * Blaschke and Cierniak (2021) D. Blaschke and M. Cierniak, Astron. Nachr. 342, 227 (2021), eprint 2012.15785.
  * Miao et al. (2020) Z. Miao, A. Li, Z. Zhu, and S. Han, Astrophys. J. 904, 103 (2020), eprint 2006.00839.
  * Paschalidis et al. (2018) V. Paschalidis, K. Yagi, D. Alvarez-Castillo, D. B. Blaschke, and A. Sedrakian, Phys. Rev. D 97, 084038 (2018), eprint 1712.00451.
  * Alford and Sedrakian (2017) M. G. Alford and A. Sedrakian, Phys. Rev. Lett. 119, 161104 (2017), eprint 1706.01592.
  * Li et al. (2020) J. J. Li, A. Sedrakian, and M. Alford, Phys. Rev. D 101, 063022 (2020), eprint 1911.00276.
  * Antić et al. (2021) S. Antić, M. Shahrbaf, D. Blaschke, and A. G. Grunfeld (2021), eprint 2105.00029.
  * Contrera et al. (2022) G. A. Contrera, D. Blaschke, J. P. Carlomagno, A. G. Grunfeld, and S. Liebing, Phys. Rev. C 105, 045808 (2022), eprint 2201.00477.
  * Antoniadis et al. (2013) J. Antoniadis et al., Science 340, 6131 (2013), eprint 1304.6875.
  * Riley et al. (2021) T. E. Riley et al., Astrophys. J. Lett. 918, L27 (2021), eprint 2105.06980.
  * Riley et al. (2019) T. E. Riley et al., Astrophys. J. Lett. 887, L21 (2019), eprint 1912.05702.
  * Raaijmakers et al. (2019) G. Raaijmakers et al., Astrophys. J. Lett. 887, L22 (2019), eprint 1912.05703.
  * Abbott et al. (2018) B. P. Abbott et al. (LIGO Scientific, Virgo), Phys. Rev. Lett. 121, 161101 (2018), eprint 1805.11581.
  * Bauswein et al. (2017) A. Bauswein, O. Just, H.-T. Janka, and N. Stergioulas, Astrophys. J. Lett. 850, L34 (2017), eprint 1710.06843.
  * Annala et al. (2018) E. Annala, T. Gorda, A. Kurkela, and A. Vuorinen, Phys. Rev. Lett. 120, 172703 (2018), eprint 1711.02644.
  * Capano et al. (2020) C. D. Capano, I. Tews, S. M. Brown, B. Margalit, S. De, S. Kumar, D. A. Brown, B. Krishnan, and S. Reddy, Nature Astron. 4, 625 (2020), eprint 1908.10352.
  * Dietrich et al. (2020) T. Dietrich, M. W. Coughlin, P. T. H. Pang, M. Bulla, J. Heinzel, L. Issa, I. Tews, and S. Antier, Science 370, 1450 (2020), eprint 2002.11355.
  * Gerlach (1968) U. H. Gerlach, Phys. Rev. 172, 1325 (1968).
  * Alvarez-Castillo et al. (2019) D. E. Alvarez-Castillo, D. B. Blaschke, A. G. Grunfeld, and V. P. Pagura, Phys. Rev. D99, 063010 (2019), eprint 1805.04105.
  * Blaschke et al. (2020a) D. Blaschke, A. Ayriyan, D. E. Alvarez-Castillo, and H. Grigorian, Universe 6, 81 (2020a), eprint 2005.02759.
  * Glendenning and Kettner (2000) N. K. Glendenning and C. Kettner, Astron. Astrophys. 353, L9 (2000), eprint astro-ph/9807155.
  * Benic et al. (2014) S. Benic, D. Blaschke, G. A. Contrera, and D. Horvatic, Phys. Rev. D 89, 016007 (2014), eprint 1306.0588.
  * Alvarez-Castillo and Blaschke (2015) D. E. Alvarez-Castillo and D. Blaschke, Phys. Part. Nucl. 46, 846 (2015), eprint 1412.8463.
  * Zacchi et al. (2017) A. Zacchi, L. Tolos, and J. Schaffner-Bielich, Phys. Rev. D 95, 103008 (2017), eprint 1612.06167.
  * Alvarez-Castillo and Blaschke (2017) D. E. Alvarez-Castillo and D. B. Blaschke, Phys. Rev. C 96, 045809 (2017), eprint 1703.02681.
  * Montana et al. (2019) G. Montana, L. Tolos, M. Hanauske, and L. Rezzolla, Phys. Rev. D 99, 103009 (2019), eprint 1811.10929.
  * Blaschke et al. (2020b) D. Blaschke, D. E. Alvarez-Castillo, A. Ayriyan, H. Grigorian, N. K. Largani, and F. Weber, in _Topics on Strong Gravity_ , edited by C. A. Vasconcellos (World Scientific, Singapore, 2020b), pp. 207–256, eprint 1906.02522.
  * Abbott et al. (2020) R. Abbott et al. (LIGO Scientific, Virgo), Astrophys. J. Lett. 896, L44 (2020), eprint 2006.12611.
  * Yudin et al. (2014) A. V. Yudin, T. L. Razinkova, D. K. Nadyozhin, and A. D. Dolgov, Astron. Lett. 40, 201 (2014), eprint 1404.0865.
  * Cierniak and Blaschke (2020) M. Cierniak and D. Blaschke, Eur. Phys. J. ST 229, 3663 (2020), eprint 2009.12353.
  * Cierniak and Blaschke (2021) M. Cierniak and D. Blaschke, Astron. Nachr. 342, 819 (2021), eprint 2106.06986.
  * Cierniak and Blaschke (2022) M. Cierniak and D. Blaschke, EPJ Web Conf. 258, 07009 (2022), eprint 2112.09166.
  * Peters and Mathews (1963) P. C. Peters and J. Mathews, Phys. Rev. 131, 435 (1963).
  * Abbott et al. (2017) B. P. Abbott et al. (LIGO Scientific, Virgo), Phys. Rev. Lett. 119, 161101 (2017), eprint 1710.05832.
  * Abbott et al. (2019) B. P. Abbott et al. (LIGO Scientific, Virgo), Phys. Rev. X 9, 031040 (2019), eprint 1811.12907.
  * Yagi and Yunes (2013) K. Yagi and N. Yunes, Phys. Rev. D 88, 023009 (2013), eprint 1303.1528.