Impact of Multiple Phase Transitions in Dense QCD on Compact Stars

Abstract

This review covers several recent developments in the physics of dense QCD with an emphasis on the impact of multiple phase transitions on astrophysical manifestations of compact stars. To motivate the multi-phase modeling of dense QCD and delineate the perspectives, we start with a discussion of the structure of its phase diagram and the arrangement of possible color-superconducting and other phases. It is conjectured that pair-correlated quark matter in $\beta$-equilibrium is within the same universality class as spin-imbalanced cold atoms and the isospin asymmetrical nucleonic matter. This then implies the emergence of phases with broken space symmetries and tri-critical (Lifshitz) points. The beyond-mean-field structure of the quark propagator and its non-trivial implications are discussed in the cases of two- and three-flavor quark matter within the Eliashberg theory, which takes into account the frequency dependence (retardation) of the gap function. We then construct an equation of state (EoS) that extends the two-phase EoS of dense quark matter within the constant speed of sound parameterization by adding a conformal fluid with a speed of sound $c_{\mathrm{conf}.}=1/\sqrt{3}$ at densities $\ge 10n_{\mathrm{sat}}$, where $n_{\mathrm{sat}}$ is the saturation density. With this input, we construct static, spherically symmetrical compact hybrid stars in the mass–radius diagram, recover such features as the twins and triplets, and show that the transition to conformal fluid leads to the spiraling-in of the tracks in this diagram. Stars on the spirals are classically unstable with respect to the radial oscillations but can be stabilized if the conversion timescale between quark and nucleonic phases at their interface is larger than the oscillation period. Finally, we review the impact of a transition from high-temperature gapped to low-temperature gapless two-flavor phase on the thermal evolution of hybrid stars.

1. Introduction

The astrophysics of compact stars entered the era of multimessenger astronomy in 2017 with the discovery of the neutron star binary merger event GW170817 [1]. Combined with radio observations of massive pulsars in binaries with white dwarfs [2] and X-ray observations of nearby solitary neutron stars [3,4], compact star astrophysics nowadays offers important insights into their global properties and potentially into the phase structure of dense matter [5,6,7,8,9,10,11,12,13]. Studies of matter at high densities are fundamental to our understanding of the strong force of the Standard Model and underlying concepts such as confinement, spontaneous chiral symmetry breaking, and dynamical mass generation [14].

This work studies the impact of multiple phase transitions in dense QCD matter on the physics of compact stars. We partially review the relevant physics but also provide a novel discussion of the mass–radius diagram of compact stars in the case where a conformal fluid is added to the two-phase, constant speed of sound parametrization of the EoS of quark matter [15]. This is motivated by the recent work that showed that even though the high-central-density stars are typically unstable toward radial oscillation (see Ref. [16], hereafter BTM), i.e., when $dM/d{\rho }_c<0$, where M is the star’s gravitational mass and ${\rho }_c$ is the central density, they can be stabilized if the conversion between nucleonic and quark phases is slow compared to the characteristic period of radial oscillations [17,18,19,20].

To motivate the modeling, Section 2 provides a brief overview of the phase diagram of dense QCD matter as we understand it from the studies of the thermodynamics of various high-density phases, such as the color-superconducting phases [6,7,8,9,10] or quarkionic phases [21,22,23,24,25]. Utilizing the knowledge gained from the studies of imbalanced cold atoms [26] and isospin asymmetrical nuclear matter [27,28], the possible phase structure of pair-correlated quark matter is conjectured based on the universal features of imbalanced pair-correlated fermionic systems. Section 3 discusses the computations of Green’s functions in the two- and three-flavor phases [29,30] and potential new effects arising beyond the adiabatic (frequency-independent) approximation of the gap. In Section 4, the constant speed-of-sound parameterization of the EoS of quark matter phases [31,32] is used to explore the mass–radius (M-R) diagram of compact stars with deconfinement and multiple phase transitions. For two-phase transitions, one from nucleonic to two-flavor quark matter and another from two-flavor to three-flavor phase of quark matter, we recover the fourth family of compact stars, which is separated from the third family by the instability region [15]. Here we show that a high-density phase of conformal fluid at densities $\ge 10n_{\mathrm{sat}}$, where $n_{\mathrm{sat}}=0.16$ fm${}^{-3}$ is the saturation density, modifies the classically unstable tracks in the M-R diagram compared to the case when such transition does not occur [33]. Such modification is phenomenologically important because of the possible stabilization mechanism of radial oscillation modes of hybrid stars [17,18,19,20], which is discussed in Section 6. In Section 5, we discuss the cooling of compact stars with quark cores. We then simulate the thermal evolution of these stars on a time scale on the order of million years with a focus on the impact of the phase transition from the gapped to the gapless phase of 2SC matter in the core of the star. Our conclusions are given in Section 7.

2. A Brief Review of the Phase Diagram of Dense QCD

Matter in compact stars covers the large number density ($n\ge n_{\mathrm{sat}}$), large isospin, and relatively low-temperature ($0\le T\le 100$ MeV) portion of the phase diagram of strongly interacting matter. The extremely low temperature ($T\le 0.1$) The MeV regime is relevant for mature compact stars, whereas the higher temperature domain is relevant for supernovas and binary neutron star mergers. The complexity of the phase diagram arises due to the multiple order parameters describing (interrelated) phenomena, which include deconfinement phase transition (with the Polyakov loop as the order parameter of the center symmetry), chiral phase transition (and its condensate as the order parameter), the color-superconducting phases (with the anomalous correlator as the order parameter). Depending on the non-zero value of one or several order parameters, distinct phases may arise: an extensively studied case is color-superconducting phases with various pairing patterns [6,7,8,9,10]. A more recent suggestion is a confined but chirally symmetric quarkyonic phase at compact star densities [21,23,24]. A crude sketch of the phase diagram of strongly interacting matter is shown Figure 1, along with the regions that are covered by current and future facilities (RHIC, NICA, and FAIR). The low-density and low-temperature region of the phase diagram contains nuclei embedded into a sea of charge-neutralizing electrons and neutrons at higher densities. As the density increases, a first-order phase transition to bulk nuclear matter occurs at around $0.5n_{\mathrm{sat}}$. A further increase in density can lead to the deconfinement of nucleons to form quark matter for $n\ge (2-3)\times n_{\mathrm{sat}}$.

The transition from nuclear matter to deconfined quark matter could be of the first or second order, or a crossover [6,7,8,9,10]. The first-order phase transition leaves a marked imprint on the macroscopic properties of compact stars because the EoS contains a density jump, which may give rise to new stable branches of compact stars (i.e., their third family) separated from nucleonic stars by a region of instability [34,35,36,37]. Smooth crossover without changes in the values of the order parameter or the wave function of the three-quark states would be a less dramatic change in the slope of the EoS, best visualized in terms of the speed of sound [10,22,38]. As mentioned above, two sequential first-order phase transitions can lead to the appearance of a new branch of compact stars—fourth family—separated from the third family by an instability region [15,33,39,40], assuming the classical stability criterion $dM/d{\rho }_c>0$ is valid. In the case of slow phase transition between the nuclear and quark matter phases, the two families are not separated by an instability region; i.e., they form a continuous branch where the regions with $dM/d{\rho }_c<0$ are stabilized [20] (see Section 6).

Deconfined quark matter at low temperature and high density is expected to have characteristic features of degenerate Fermi systems, which are familiar from condensed matter physics. Therefore, emergent phenomena such as (color) superconductivity are expected in channels where gluon exchange is attractive [6,7,8]. Various color superconducting phases may arise, depending on the number of flavors and colors involved in the pairing, the ratio of $\Delta /\delta \mu$, where $\Delta$ is the gap in the quasiparticle spectrum, $\delta \mu =({\mu }_d-{\mu }_d)/2$ is the difference in the chemical potentials of down (d) and up (u) quarks, and the mass of strange quark $m_s$ in the three-flavor quark matter. The two-flavor candidate phases classified according to these parameters are as follows:

(a)

The “2SC” phase (where the abbreviation refers to two-superconducting-colors) [6]

$$\begin{array}{c}\hfill {\Delta }_{2SC}\propto \langle {\psi }^T\left(x\right)C{\gamma }_5{\tau }_2{\lambda }_2\psi \left(x\right)\rangle \ne 0,0\le \delta \mu <\Delta /\sqrt{2},\end{array}$$

(1)

where $C=i{\gamma }^2{\gamma }^0$ is the charge conjugation operator, ${\tau }_2$ is the second component of the Pauli matrix acting in the SU${\left(2\right)}_f$ flavor space, and ${\lambda }_A$ is the antisymmetric Gell–Mann matrix acting in the SU${\left(3\right)}_c$ color space. The properties of the 2SC phase resemble that of the ordinary BCS theory, including vanishing resistivity and vanishing heat capacity because the quarks near the Fermi surface remain gapped.

(b)

Phases with broken space symmetries, which are associated with a finite momentum of the condensate [6,8] (hereafter FF phase) or deformation of the Fermi surface [41] (hereafter DFS phase):

$$\begin{array}{ccc}\hfill & & {\Delta }_{2SC}\ne 0,\delta \mu >\Delta /\sqrt{2},\overrightarrow{P}\ne 0\left(\mathrm{FF}\right)\hfill \end{array}$$

(2)

$$\begin{array}{ccc}& & {\Delta }_{2SC}\ne 0,\delta \mu >\Delta /\sqrt{2},\delta ϵ\ne 0\left(\mathrm{DFS}\right)\hfill \end{array}$$

(3)

where $\overrightarrow{P}$ is the center of mass momentum of a Cooper pair and $\delta ϵ$ quantifies the quadrupole deformation of the Fermi surfaces of u and d quarks.

(c)

Mixed phase(s) [42]

$$\begin{array}{c}\hfill {\Delta }_{2SC}\propto \langle {\psi }^T\left(x\right)C{\gamma }_5{\tau }_2{\lambda }_2\psi \left(x\right)\rangle \ne 0,\delta \mu =0,0\le x_s\le 1,\end{array}$$

(4)

which corresponds to a mixture between a perfectly symmetrical “2SC” superconductor and a normal system accommodating the excess number of d quarks. Here, $x_s$ is the filling factor defined as the ratio of the superconducting and total volumes. (We assume that there is an excess of d over u quarks, as is expected in quark matter in compact stars under $\beta$-equilibrium.)

The color-flavor-locked (CFL) phase [43] is expected to be the ground state of three-flavor quark matter at asymptotically large densities where the strange quark is massless, Fermi surfaces of quarks coincide, and, therefore, the pairing among quarks occurs in a particularly symmetrical manner. At densities relevant for neutron stars, the perfect CFL phase is unlikely to be realized; rather, some of its variants have $m_s\ne 0$ and/or $\delta \mu \ne 0$—chemical potential shifts between various flavors of quarks [44]. Therefore, the phases listed above can be replicated with an allowance of additional non-zero $us$ and $ds$ pairings

$$\begin{array}{c}\hfill {\Delta }_{ud}\ne 0,{\Delta }_{sd}\ne 0,{\Delta }_{su}\ne 0,(m_s\ne 0;\delta \mu \ne 0).\end{array}$$

(5)

A complete phase diagram of quark matter that includes most, if not all, of the phases mentioned above, is not available to date. However, various imbalanced superfluids, such as cold atoms, isospin asymmetrical nuclear matter, and flavor-imbalanced quark matter show a high degree of universality. Thus, possible structures of the phase diagram of quark matter can be conjectured by extrapolating from the detailed studies of the phase diagrams of cold atomic gases [26] and isospin asymmetrical nuclear matter [27]. These are, clearly, speculative and need to be confirmed using explicit computations of relevant quark phases.

Figure 2 shows two schematic phase diagrams of color-superconducting matter in the density–temperature plane. For sufficiently large temperatures, the unpaired normal phase is the preferred state of matter, ignoring any other correlation beyond the pairing. The phases with broken symmetries, the FF and the DFS phases, are preferable in temperature–density strips at low temperatures and high densities. At lower temperatures, the PS phase is the preferred one. At higher temperatures, the spatially symmetric 2SC phase dominates. It is seen that the phase diagram contains two tri-critical points, i.e., the points where three different phases coexist. The critical point, which has the FF state at the intersection, is a Lifshitz point as, per construction, it is a meeting point of the modulated (FF), ordered (PS/2SC), and disordered (unpaired) states. Of course, this is the case if the transition temperature to the CFL phase is below the tri-critical temperature; otherwise, the unpaired state should be replaced by a variant of the CFL phase. Note that depending on the parameters of the model, two or one tricritical points may be located on the unpairing line or the line of transition to the CFL phase, as illustrated in Figure 2, left and right panels, respectively. The model can be tuned to produce a four-critical point if both points coincide. We also note that the low-density limit corresponds to the strong coupling regime where the pairs are tightly bound, whereas the high-density limit corresponds to the weak-coupling regime. Therefore, one can anticipate signatures of BCS–BEC crossover. These can be seen by examining several characteristic quantities, for example, the ratio of the coherence length to the interparticle distance $\xi /d$, where $\xi /d\gg 1$ corresponds to the BCS and $\xi /d\ll 1$ corresponds to the BEC limit, or the ratio of the gap to the (average) chemical potential $\Delta /\mu$, where $\Delta /\mu \ll 1$ corresponds to the BCS and $\Delta /\mu \gg 1$ corresponds to the BEC limit. For discussions of BCS=-BEC crossover in dense quark matter, see Refs. [45,46,47,48]. This phenomenon shows a high degree of universality as well; see for example, the studies of nuclear matter [49,50,51] and cold atoms [26,52].

3. Structure of Green Functions in Two-Flavor Quark Matter

Order-by-order computations of the magnitude of the gaps in the superconducting phases can be carried out in the weak-coupling (extreme high-density) regime, where the one-gluon exchange is the dominant interaction. Approximate Eliashberg-type equations for the flavor-symmetric 2SC phase were solved within one-gluon exchange approximation in Refs. [53,54], showing that the pairing gap scales with the coupling g as $\Delta \sim \mu g^{-5}exp(-1/g)$. Such a scaling also applies to the high-density CFL phase, where the perturbative approach is more reliable than at densities relevant to the 2SC phase. More recently, Eliashberg-type equations were solved for two-flavor [29] and three-flavor superconductors [30]. The first study used the quark–meson coupling model, keeping only the frequency dependence of the gap, whereas the second study kept frequency and momentum dependences but ignored the imaginary part of the pairing gap. These theories not only improve the description of quark matter but also lead to phenomenologically important implications, such as the presence of electrons in the CFL phase [30], which are not allowed when the gap is constant [55].

We now briefly outline these approaches following Ref. [29]. The inverse Nambu–Gorkov quark propagator is given by

$$S^{-1}\left(q\right)=\left(\begin{array}{cc}\cancel{q}+\mu {\gamma }_0-m& \overline{\Delta }\\ \Delta & {(\cancel{q}-\mu {\gamma }_0+m)}^T\end{array}\right),$$

(6)

where q is the four-momentum and $\Delta$ is the gap with $\overline{\Delta }\equiv {\gamma }_0{\Delta }^{†}{\gamma }_0$. Equation (6) is written for the case of equal number densities of up and down quarks with a common chemical potential $\mu$ and mass m. The bare quark–meson vertices ${\Gamma }_{\pi }^i\left(q\right)$ and ${\Gamma }_{\sigma }\left(q\right)$ are given by

$$\begin{array}{c}\hfill {\Gamma }_{\pi }^i\left(q\right)=\left(\begin{array}{cc}\frac{{\tau }^i}{2}{\gamma }_5& 0\\ 0& -{\left(\frac{{\tau }^i}{2}{\gamma }_5\right)}^T\end{array}\right),{\Gamma }_{\sigma }\left(q\right)=\left(\begin{array}{cc}\mathbb{I}& 0\\ 0& -\mathbb{I}\end{array}\right),\end{array}$$

(7)

where pions couple to quarks using a pseudo-scalar coupling, whereas $\sigma$s couple via a scalar coupling, with $\mathbb{I}$ being a unit matrix in the Dirac and isospin spaces. Their propagators are given by

$$\begin{array}{c}\hfill D_{\pi }\left(q\right)=\frac{1}{q_0^2-{\mathit{q}}^2-m_{\pi }^2},D_{\sigma }\left(q\right)=\frac{1}{q_0^2-{\mathit{q}}^2-m_{\sigma }^2},\end{array}$$

(8)

where $m_{\pi /\sigma }$ values are the meson masses. The equation for the gap in the Fock approximation is given via

$$\begin{array}{ccc}\hfill \Delta \left(k\right)& =& i{g}_{\pi }^2\int \frac{d^{4}q}{{\left(2\pi \right)}^4}{\left(-\frac{{\tau }^i}{2}{\gamma }_5\right)}^T{S}_{21}\left(q\right)\frac{{\tau }^j}{2}{\gamma }_5{\delta }_{ij}{D}_{\pi }(q-k)\hfill \\ & +& i{g}_{\sigma }^2\int \frac{d^{4}q}{{\left(2\pi \right)}^4}{(-\mathbb{I})}^T{S}_{21}\left(q\right)\mathbb{I}{D}_{\sigma }(q-k),\hfill \end{array}$$

(9)

where $g_{\pi }$ and $g_{\sigma }$ are the coupling constants. Adopting the color–flavor structure of the gap function corresponding to a 2SC superconductor, one then finds

$$\begin{array}{c}\hfill {\Delta }_{ij}^{ab}\left(k\right)={\left({\lambda }_2\right)}^{ab}{\left({\tau }_2\right)}_{ij}C{\gamma }_5[{\Delta }_{+}\left(k\right){\Lambda }^{+}\left(k\right)+{\Delta }_{-}\left(k\right){\Lambda }^{-}\left(k\right)],\end{array}$$

(10)

where $a,b\cdots$ refer to the color space, $i,j,\cdots$ refer to the flavor space, and the projectors onto the positive and negative states are defined in the standard fashion as ${\Lambda }^{\pm }\left(k\right)=(E_k^{\pm }+\mathit{\alpha }·\mathit{k}+m{\gamma }_0)/2E_k^{\pm },$ where $E_k^{\pm }=\pm \sqrt{{\mathit{k}}^2+m^2}$ and $\mathit{\alpha }={\gamma }_0\mathit{\gamma }$. The coupled Equations (6)–(10) must be solved for the gap function, which is a function of three-momentum and the frequency. In the low-temperature limit, the relevant momenta are close to the Fermi momentum and the dependence on the magnitude of the three-momentum can be eliminated by fixing it at the Fermi momentum. The gap Equation (9) then depends only on the energy, which reflects the fact that the pairing interaction is not instantaneous—a common feature of the Fock self-energies in ordinary many-body perturbation theory. The solutions for the positive energy projection of the gap function are shown in Figure 3 as a function of frequency. The structure of the real and imaginary components of the gap function shows a maximum around frequencies at which the meson spectral functions are peaked. Thus, it is important to include the retardation effect when the color superconductor is probed at such frequencies. In the low-frequency limit, it is sufficient to use the BCS approximation where the interaction is instantaneous so that the imaginary part vanishes Im$\Delta \left(\omega \right)=0$ and the real part is a constant Re$\Delta =\Delta (\omega =0)$.

Ref. [30] considered the full momentum and energy dependence of the gap in the Fock approximation within the Yukawa model but neglected the imaginary part of the anomalous self-energy. Their work shows that the retardation implies a CFL phase that is not a perfect insulator, as charge neutrality requires some electrons to be present in matter. This is not the case in the treatment based on the BCS model [55]. Thus, the phenomenology of the CFL phase is modified: its specific heat, thermal conductivity, and magnetic response will change due to the contribution of electrons. This example, in which the simple BCS ansatz for the gap is replaced by a more complete gap function, demonstrates some unexpected features of color superconductors, which may be important for their transport and dynamic response to various probes.

4. Equation of State and Mass–Radius Diagram

We have seen in the previous section that the phase diagram of quark matter may have a complicated structure. At minimum, there are two robust phases of color superconducting matter: the low-density two-flavor 2SC phase and the high-density three-flavor CFL phase. See, for example, Ref. [56] for a Nambu–Jona-Lasinio study and compact stars with two phase transitions in this model. However, additional phases are very likely because it is energetically favorable to break the rotational and translational symmetries due to the stress on the paired state induced by the finite mass of the strange quark and $\beta$-equilibrium, which induces disparity in the chemical potentials of u and d quarks. In addition, or alternatively, quarkyonic phases may interfere.

For the specific computation below, we adopt a covariant density functional EoS of nuclear matter in the nucleonic phase [39,40,57]. This EoS, in the absence of the phase transition to quark matter, produces nucleonic compact stars with a maximum $m\equiv M/M_{\odot }\simeq 2.6$, where M is the gravitational mass of the star and $M_{\odot }$ is the solar mass. Allowing for phase transition to quark matter, we consider a straightforward extension of the constant speed of sound EoS of Ref. [15] that allows for a conformal phase of quark matter at high densities with a constant speed of sound, i.e.,

$$\begin{array}{c}\hfill p\left(ϵ\right)=\left\{\begin{array}{cc}{p}_1,\hfill & {ϵ}_1<ϵ<{ϵ}_1+\Delta {ϵ}_1\hfill \\ p_1+s_1\left[ϵ-({ϵ}_1+\Delta {ϵ}_1)\right],\hfill & {ϵ}_1+\Delta {ϵ}_1<ϵ<{ϵ}_2\hfill \\ p_2,\hfill & {ϵ}_2<ϵ<{ϵ}_2+\Delta {ϵ}_2\hfill \\ p_2+s_2\left[ϵ-({ϵ}_2+\Delta {ϵ}_2)\right],\hfill & {ϵ}_2+\Delta {ϵ}_2<ϵ<{ϵ}_3\hfill \\ p_3,\hfill & {ϵ}_3<ϵ<{ϵ}_3+\Delta {ϵ}_3\hfill \\ p_3+s_3\left[ϵ-({ϵ}_3+\Delta {ϵ}_3)\right],\hfill & ϵ>{ϵ}_3\hfill \end{array}\right.\end{array}$$

(11)

where the three pairs of the pressure and energy density $p_{1,2,3}$ and ${ϵ}_{1,2,3}$ correspond, respectively, to the transition from hadronic to quark matter, from a low-density (2SC) quark phase to a high-density (CFL) quark phase, and from the high-density quark phase to the conformal fluid. The squared sound speeds in the quark phases are denoted by $s_1$, $s_2$, and $s_3=c_{\mathrm{conf}.}^2$. Note that we assume that the 2SC and CFL quark phases are separated by a jump at the phase boundary, as it follows from the study of Ref. [56]. At high densities, the CFL phase reaches the “conformal limit” where the interactions are dominated by the underlying conformal symmetry of QCD. In this limit, the speed of sound squared is $s_3=1/3$ (in units of speed of light), whereas the effects of the pairing gap of the CFL phase can be neglected in a first approximation. Note that we allow for a small jump between proper CFL and conformal zero-gap fluid, but its effect on the observables is marginal, i.e., a smooth interpolation would not change the results.

According to Equation (11), the modeling of the EoS of quark phases involves the following parameters:

• The three (energy) densities at which the sequential transitions between the nucleonic phase, 2SC phase, CFL phase, and conformal fluid take place.
• The magnitudes of the jumps in the energy density at the points of the transition from nuclear to the 2SC phase, $\Delta {\epsilon }_1$, from the 2SC to the CFL phase, $\Delta {\epsilon }_2$, and from the CFL to the conformal fluid phase $\Delta {\epsilon }_3$.
• The speeds of sound in the 2SC and CFL phases $s_1$ and $s_2$. The speed of sound of the conformal fluid is held fixed at $s_3=1/3$. Note that for any phase, $s\le 1$ by causality.

Our model EoS is constructed using the following parameters. The transition pressure and energy density from nuclear and quark matter are $p_1=1.7\times {10}^{35}$ dyn/cm${}^2$ and ${ϵ}^{}1=8.4\times {10}^{14}$ g cm${}^{-3}$, respectively. The magnitude of the first jump $\Delta {ϵ}^{}1=0.6{ϵ}^{}1$. The upper range of the energy density of the 2SC phase is determined as ${ϵ}^{}{1}^{\max }={\delta }_{2SC}({ϵ}^{}1+\Delta {ϵ}^{}1)$, where ${\delta }_{2SC}$ is a dimensionless parameter measuring the width of the 2SC phase. The magnitude of the second jump is parametrized in terms of the ratio parameter $r=\Delta {ϵ}^{}2/\Delta {ϵ}^{}1$. The extent of the CFL phase is determined by limiting its energy density range to ${ϵ}^{}{2}^{\max }={\delta }_{CFL}({ϵ}^{}2+\Delta {ϵ}^{}2)$. The transition to the conformal fluid is assumed to be of the first order with a small (compared to other scales) energy-density jump equal to $0.1r$. The transition to the conformal fluid phase occurs at densities ${ϵ}^{}3=2.25$–$2.57\times {10}^{15}$ g cm${}^{-3}$, i.e., by about a factor of 10 larger than the saturation density. The speeds of sound squared are fixed as

$$\begin{array}{c}\hfill s_1=0.7,s_2=1,s_3=\frac{1}{3}.\end{array}$$

(12)

The values of $s_1$ and $s_2$ are chosen to obtain triplet configurations with large enough masses of hybrid stars. The magnitudes of jumps between the nuclear, 2SC, and CFL phases were chosen suitably to produce twin and triplet configuations [15].

Figure 4 shows a collection of EoS constructed based on Equation (11), which shares the same low-density nuclear EoS. In this collection, we vary the parameter r (as indicated in the plot) for fixed values ${\delta }_{2SC}={\delta }_{CFL}=0.27$. The corresponding M–R relations for static, spherically symmetrical stars obtained from solutions to Tolman–Oppenheimer–Volkoff equations are shown in Figure 5.

For the chosen magnitude of the first jump $\Delta {ϵ}^{}1$, the M–R curves show the phenomenon of twins—two stars of the same mass but different radii. The radii of twins differ by about 1 km. The more compact configuration is a hybrid star, i.e., a star with a quark core and nuclear envelope, whereas the less compact counterpart is a purely nucleonic star. The second phase transition may or may not result in a classically stable sequence depending on the value of the parameter r parameterizing the magnitude of the second jump. For small jumps $r=0.1$ and 0.23, new stable branches arise, which are continuously connecting to the stable 2SC branch ($r=0.1$) or are separated by a region where the stars are unstable ($r=0.23$). It is seen that, in this case, triplets of stars with different radii but the same masses appear. The densest stars contain, in addition to the 2SC phase, a layer of the CFL phase, whereby the central density on the stable branch can exceed the onset density of the conformal fluid. This implies that the densest member of a triplet will contain in its center conformal fluid with $c_{\mathrm{conf}.}=1/\sqrt{3}$. For each M-R curve in Figure 5, the star with a central density at which the conformal fluid first appears is shown by a dot (this density is fixed at $10n_{\mathrm{sat}}$). The stable branch of conformal fluid containing stars is followed by a classically unusable branch with $dM/d{\rho }_c<0$. For asymptotically large central densities, the masses and radii increase again. The family of the EoSs that differ only in the value of the parameter r cross at a “special point”. This type of crossing has been observed for twin star configurations with a variation in a particular parameter of the EoS [58]; however, the EoS excluded two sequential phase transitions. The behavior of M–R curves at very high central densities differs from the ones that were found in Ref. [33], where a branch of ultracompact twin stars with masses of the order of $1M_{\odot }$ and radii in the range of 6–7 km were found for a single phase transition from the nuclear matter to the quark phase. Thus, we conclude that the high-density asymptotics of the EoS modifies the behavior of the M–R curves if the conformal limit is achieved at densities of the order of $10n_{\mathrm{sat}}$.

The observation above may have phenomenological implications for the following reason. The stability of stellar configurations is commonly determined by the requirement that the star’s mass must increase with increasing central density (or central pressure), i.e., $\partial M/\partial {\rho }_c>0$. An alternative and physically more transparent method is to compute the radial modes of oscillation of a star and determine the stable configurations from the requirement that their frequencies are real. Ref. [17] showed that the classical stability conditions fail if the conversion rate is slow, i.e., if its characteristic timescale is longer than the period of oscillations. In that case, the fundamental modes are stable even when $\partial M/\partial {\rho }_c<0$; i.e., stars with central densities larger than the one corresponding to the maximum-mass star (which lie to the left from the maximum on the M–R diagram in Figure 5) will be stable. This observation also applies to configurations with two-phase transitions, as shown in Refs. [19,20]. Furthermore, Ref. [20] shows that, in this case, the classically unstable stars contribute to the count of same-mass stars, which leads to the appearance of higher-order multiplets such as quadruplets, quintuplets, and sextuplets. We will return to the stability of hybrid stars in Section 6.

5. Cooling of Compact Stars with Quark Matter Cores

The cooling of compact stars may provide indirect information about quark phases in hybrid stars. The properties of phases of dense quark matter affect both neutrino emission and the specific heat content that determine the cooling rate of a compact object in general; see Refs. [59,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75].

Non-superconducting relativistic quark matter cools predominantly via the direct Urca processes involving d, u, and s quarks [76]

$$\begin{array}{c}\hfill \begin{array}{cc}& d\to u+e+{\overline{\nu }}_e,\hfill \\ & u+e\to d+{\nu }_e,\hfill \\ & s\to u+e+{\overline{\nu }}_e,\hfill \\ & u+e\to s+{\nu }_e,\hfill \end{array}\end{array}$$

(13)

where ${\nu }_e$ and ${\overline{\nu }}_e$ are the electon neutrino and antineutrions. The neutrino emissivity through the direct Urca process for non-strange quarks is given via [77]

$$\begin{array}{c}\hfill {ϵ}_{\beta }=8.8\times {10}^{26}{\alpha }_s\left(\frac{n}{n_{\mathrm{sat}}}\right)Y_e^{1/3}{T}_9^6\mathrm{ergs}{\mathrm{cm}}^{-3}{\mathrm{s}}^{-1},\end{array}$$

(14)

where n is the baryon density, $Y_e$ is the electron fraction, $T_9$ is the temperature in units of ${10}^9\mathrm{K}$, and ${\alpha }_s$ is the running strong coupling constant. The emissivity given by Equation (14) implies that the stars containing unpaired quark matter would cool quickly via this direct Urca process. The cooling would be slower if the quark spectrum contains a gap. In the case of the phenomenologically relevant 2SC phase, two alternatives are possible, depending on whether the Fermi surfaces of quarks are a full gap or they contain zero-gap segments (nodes). The latter feature arises in the case of pairing between fermions on different Fermi surfaces, as discussed in Section 2.

Ref. [78] studied a generic case where the quark spectrum is gapped if the parameter $\zeta ={\Delta }_0/\delta \mu$ associated with the new scale $\delta \mu =\left({\mu }_d-{\mu }_u\right)/2$, where ${\mu }_{u,d}$ are the chemical potentials of light quarks and ${\Delta }_0$ is the gap for $\delta \mu =0$. The suppression of emissivity by pairing is qualitatively different in the cases in which $\zeta >1$ and $\zeta <1$. The novelty arises in the second case, where Fermi surfaces have nodes and particles can be excited around these nodes without any energy cost (which is not the case for gapped Fermi surfaces). Note that in the case of the FF phase, the shift in the chemical potential is replaced by a more general function—the anti-symmetric in the flavor part of the single particle spectrum of up and down quarks. This new physics can be captured by adopting a generic parameterization of the suppression factor of the quark Urca process with pairing suggested in Ref. [78]. The neutrino emissivity of the 2SC phase ${ϵ}_{2SC}^{rg}$ can be related to the Urca rate in the normal phase (14) as

$$\begin{array}{c}\hfill {ϵ}_{2SC}^{rg}\left(\zeta ;T\le T_c\right)=2f\left(\zeta \right){ϵ}_{\beta },f\left(\zeta \right)=\frac{1}{exp\left[(\zeta -1)\frac{\delta \mu }{T}-1\right]+1},\end{array}$$

(15)

where the parameters $\zeta$ and $\delta \mu$ were introduced above, T is the temperature, and $T_c$ is the critical temperature of the phase transition from normal to the 2SC phase. Furthermore, the parameter $\zeta \left(T\right)$ is temperature-dependent and we adopt the parametrization

$$\begin{array}{c}\hfill \zeta \left(T\right)={\zeta }_i-\Delta \zeta g\left(T\right),\end{array}$$

(16)

where ${\zeta }_i$ is the initial value, $\Delta \zeta$ is the constant change in this function, and the function $g\left(T\right)$ describes the transition from the initial value ${\zeta }_i$ to the asymptotic final value ${\zeta }_f={\zeta }_i-\Delta \zeta$. The transition is conveniently modeled by the following function

$$\begin{array}{c}\hfill g\left(T\right)=\frac{1}{exp\left(\frac{T-T*}{w}\right)+1},\end{array}$$

(17)

which allows one to control the temperature of transition by adjusting the parameter $T*$ and the smoothness of the transition via the width parameter w. An additional issue to address is the role of the blue quarks that do not participate in the 2SC pairing. Blue quarks may pair among themselves due to the attractive component of the strong force as in the ordinary BCS case (as both members of the Cooper pair are on the same Fermi surface). Then, the emissivity of blue quarks in the superfluid state is given by

$$\begin{array}{c}\hfill {ϵ}_{BCS}^b\left(T\le T_{cb}\right)\simeq {ϵ}_{\beta }^b\left(T>T_{cb}\right)exp\left(-\frac{{\Delta }_b}{T}\right),\end{array}$$

(18)

where ${\Delta }_b$ is the gap in the blue quark spectrum, $T_{cb}$ is the corresponding critical temperature, and ${ϵ}_{\beta }^b$ is the neutrino emissivity of blue quarks in the normal state. As discussed in Section 4, the densest members of the triplets contain cores of CFL matter that is fully gapped. In this case, the excitations are the Goldstone modes of the CFL phase. Their emissivity, as well as the specific heat, is rather small compared to other phases due to their very small number density [79]. In the following discussion, we will ignore the role of the CFL phase in the cooling of hybrid stars. In the conformal fluid phase, we expect three-flavor pairing gap $\Delta \sim \mu g^{-5}exp(-1/g)$, $g=\sqrt{4\pi {\alpha }_s}$, with a spin–flavor structure of the CFL phase.

Let us turn to the cooling simulations of hybrid stars with a gapless 2SC superconductor. The cooling tracks are shown in Figure 6, and the input physics beyond the emissivities is discussed elsewhere [61,62,67,68]. The key parameter regulating the behavior of the cooling curves in Figure 6 is the temperature $T*$, which controls the transition from the gapped to ungapped 2SC phase. Similar results were obtained in the context of rapid cooling of the compact star in Cassiopeia (Cas) A remnant in Refs. [61,62,68]. The model has a second parameter, the gap for blue-colored quarks ${\Delta }_b$, which prohibits rapid cooling via the Urca process involving only blue quarks. The third parameter w in Equation (17) accounts for the finite time scale of the phase transition—see Refs. [62,68]—but it is important only for the fine-tuning of the cooling curves close to the age of the Cas A. The various cooling tracks shown in Figure 6 correspond to various values of $T*$ for fixed values of w and ${\Delta }_b$ and stellar configuration of mass $1.93M_{\odot }$. It is seen that if $T*$ is small, then the quark core does not influence the cooling, because during the entire evolution $T>T*$; therefore the neutrino emission is suppressed by the fully gapped Fermi surfaces of red-green quarks. For large $T*$, early transition to the gapless phase occurs, and the star cools fast via the direct Urca process. Note that the value of $T*$ can be fine-tuned to reproduce not only the current temperature of Cas A but also the fast decline claimed to be observed during the last decade or so; see Ref. [80] and references therein. From the brief discussion above, one may conclude that the phase transitions within the cold QCD phase diagram may induce interesting and phenomenologically relevant changes in the cooling behavior of compact stars. Although we will not discuss in any depth the dependence of cooling tracks on the stellar mass, it should be pointed out that the onset of new phases in the interiors of compact stars, for example, hyperonization, meson condensation, and phase transition to quark matter, lead to mass hierarchy in the cooling curves [81,82,83,84]. Typically, one finds that heavier stars that have central densities beyond the threshold for the onset of the new phase cool faster than the light stars containing only nucleonic degrees of freedom. This is also the case for models of stars studied here. For example, stars with masses $M\sim 1.1–1.6M_{\odot }$ remain warm over longer time scales and are thus hotter than their heavy analogs, which develop large quark cores.

6. Stability Criteria for Hybrid Stars

The oscillation modes of a compact star are important probes of their internal structure, as has been shown in the case of g modes, which are sensitive to the size of the density jumps at a first-order phase transition between hadronic and quark matter [85,86,87]. They are expected to leave an imprint on the emitted gravitational wave signal during the binary inspiral of a neutron star, as well as in the post-merger phase [88,89,90,91].

As discussed briefly in Section 4, high-central density stars on the descending branch of M–R diagram can have phenomenological implications if they are stabilized by some mechanism, which we discuss in this section. The main mode of instability for non-rotating, spherically symmetrical fluid stars in general relativity is the instability against the radial f-mode of oscillations [92]. If the f-mode frequency ${\omega }_f^2>0$, the stellar configuration is dynamically stable, and it is unstable if ${\omega }_f^2<0$. The location of this instability point on the M–R diagram agrees well with the turning point of the mass–central-density ($M-{\rho }_c$) curve. The stars on the ascending branch are stable, whereas those on the descending branch are unusable. The maximum mass is the point of marginal stability. Numerical simulations found some violations of this criterion [93,94], but quantitative deviations are insignificant. However, recent work found that the agreement between these criteria is strongly violated for stars with first-order phase transition, as we review below.

Early work on stellar oscillations with phase transitions inside the star was carried out in the Newtonian theory assuming uniform phases [95,96]. Two possibilities arise depending on the interplay between the scales in the problem: (a) when the conversion rate from one phase to another is fast, the interface between phases oscillates as a whole when perturbed; (b) if, however, conversion is slow, then the interface is fixed over the period of characteristic oscillations. The second case is interesting because, as shown in Ref. [17], the sign of ${\omega }_f^2$ does not change at the maximum mass $M_{\max }$ but stays positive over a segment where $\partial M/\partial {\rho }_c<0$. This implies that the classically unstable branch becomes stable against f-mode oscillations. Several subsequent studies confirmed this feature in the case of single- [18] and two-phase transitions [19,20]. The case of two-phase transition was extended in several directions in Ref. [20] by focusing on EoS, which supported classical twin and triplet star configurations, as discussed in Section 2. It was shown that in the case of slow conversion, higher-order multiplet stars arise, since now the stars on the $\partial M/\partial {\rho }_c<0$ segments of the mass–central-density curve are located on the stable branch. Also, the properties of the reaction mode of a compact star [96], which arises in case (a) with one or more rapid phase transitions, were studied.

The fundamental modes of hybrid stars are obtained from the set of equations [97,98]

$$\begin{array}{ccc}\hfill \frac{\mathrm{d}\xi }{\mathrm{d}r}& =& \left(\frac{\mathrm{d}\nu }{\mathrm{d}r}-\frac{3}{r}\right)\xi -\frac{\Delta P}{r\Gamma P},\hfill \\ \hfill \frac{\mathrm{d}\Delta P}{\mathrm{d}r}& =& \left[{\mathrm{e}}^{2\lambda }\left({\omega }^2{\mathrm{e}}^{-2\nu }-8\pi P\right)+\frac{\mathrm{d}\nu }{\mathrm{d}r}\left(\frac{4}{r}+\frac{\mathrm{d}\nu }{\mathrm{d}r}\right)\right]\left(\rho +P\right)r\xi \hfill \end{array}$$

(19)

$$\begin{array}{ccc}& -& \left[\frac{\mathrm{d}\nu }{\mathrm{d}r}+4\pi (\rho +P)r{\mathrm{e}}^{2\lambda }\right]\Delta P,\hfill \end{array}$$

(20)

where $\xi ={\xi }_{\dim }/r$, with ${\xi }_{\dim }$ being the Lagrangian displacement, r the radial coordinate, $\Delta P$ the Lagrangian perturbation of pressure, $\rho$ the mass–energy-density, $\omega$ the angular frequency, $\Gamma$ the adiabatic index, and ${\mathrm{e}}^{2\nu }$ and ${\mathrm{e}}^{2\lambda }$ the metric coefficients entering the Tomann–Oppenheimer–Volkoff equations. In a first approximation, the adiabatic index for a chemically equilibrated relativistic fluid can be taken as that of the matter in $\beta$-equilibrium $\Gamma =[(\rho +P)/P]\left(dP/d\rho \right).$ The set of Equations (19) and (20) can be solved provided the boundary conditions are known. These are specified by assuming that the displacement field is divergence-free at the center and that the Lagrangian variation of the pressure vanishes at the surface of the star:

$$\Delta P(r=0)=-3\Gamma P\xi (r=0),\Delta P(r=R)=0.$$

(21)

The ${\omega }^2$ values obtained in this manner are usually labeled according to the number of radial nodes in $\xi$ and the f mode corresponds to the nodeless mode.

In the case of multiple phase transitions in the QCD phase diagram, one needs junction conditions that relate the values of Lagrangian perturbations on both sides of the interface between phases. Such junction conditions already appear in the work of Ref. [96] in the Newtonian cases, whereas the the general relativistic case is treated in Ref. [99]. For the slow conversion rate one has the junction condition

$${\left[\Delta P\right]}_{-}^{+}=0,{\left[\xi \right]}_{-}^{+}=0;$$

(22)

for rapid conversion rate, one has

$${\left[\Delta P\right]}_{-}^{+}=0,{\left[\xi -\frac{\Delta P}{r}{\left(\frac{\mathrm{d}P}{\mathrm{d}r}\right)}^{-1}\right]}_{-}^{+}=0,$$

(23)

where $+/-$ refer to the high- and low-density sides of the transition, respectively. At present, it is not possible to state with confidence which limit is realized in quark matter, as the conversion rate varies significantly over the parameter space; see Ref. [100] for a discussion and earlier references. Ref. [20] considered modified junction conditions that smoothly interpolate between the two limiting cases.

Phenomenologically, the most interesting implication of the modified stability criteria is the existence of new stable configurations beyond those that are classically stable. In particular, in the case where twins and triplets exist according to classical criteria of stability, additional configurations will arise when conversion between phases at the interface is slow. These can form quadruplets (the maximum number in the case of twins) and quintuplets and sextuplets in the case of triplets. A particular case that allows for classical triplet stars is illustrated in Figure 7, adapted from Ref. [20]. The fundamental mode frequency ${\omega }_f$ is shown as a function of the central pressure of the stars in two cases when both interfaces (i.e., nucleonic to 2SC and 2SC-CFL) feature rapid or slow conversion. (The case of rapid–slow and slow–rapid conversions are intermediate cases, and we omitted them.) To recover the classical case, one needs to assume that the conversion at each interface is rapid: in this case, the instability region is characterized by the vanishing of the real part of ${\omega }_f$, as seen in Figure 7. In the case of slow conversions at both interfaces, one finds a continuous positive solution across the values of central densities of the stellar sequences, thus indicating that the stars are always stable, even on the descending branch of the mass–central-pressure curve.

To summarize the recent findings regarding the stability of the hybrid stars, we have seen that their stability against the fundamental oscillation modes strongly depends on the junction conditions at the interfaces between the phases. These are determined by the rate of conversion between phases at the phase boundary. In the case of slow phase transitions (i.e., when the conversion time scale is larger than the characteristic period of the oscillations), the usual stability criteria are modified and new stable segments appear that were previously unstable. Alternative variants of junction conditions that are intermediate between slow and rapid conversion were also considered, but the resulting radial modes do not differ significantly from the slow conversion case, with corresponding implications for the stability of the stars [20].

7. Conclusions

The investigation of dense QCD through the astrophysics of compact stars is an actively pursued subject. This is due to the substantial observational progress, which includes measurements of the masses and radii of pulsars and gravitational wave signals from mergers of two neutron stars and neutron-star–black-hole binaries. A more thorough comprehension of the thermodynamics of dense QCD, weak interactions, and the dynamics of phase transitions would greatly enhance our ability to model astrophysical phenomena relevant to current observational programs.

This work gave an overview of the phase diagram of cold and dense QCD appropriate for compact stars. We stressed that the universality of the phase diagram of imbalanced fermionic superfluids, such as cold atomic gases and nuclear matter, provides a valuable guide to the possible arrangement of the color-superconducting phases in neutron stars, the presence of tri-critical points, and BCS–BEC crossovers. The universality allows one to conjecture the possible structures of the phase diagram in the density–temperature plane including such phases, such as the Fulde–Ferrel phase, deformed Fermi surface phase, and the phase separation.

As a novel contribution, the previously proposed parametrization of the EoS of dense quark matter with sequential phase transitions was extended to include a conformal fluid at large densities ($n\ge 10n_{\mathrm{sat}}$) with the speed of sound $c_{\mathrm{conf}.}=1/\sqrt{3}$. The part of the M-R diagram that contains twins and triplets remains intact because the transition to conformal fluid occurs at larger central densities than those achieved in these objects. Nevertheless, for large central densities, we find behavior that is qualitatively different from earlier studies of this regime: the M–R curves spiral in; i.e., after reaching a minimum, they turn to the right (larger radius region), thus avoiding the region of ultra-compact stars. Therefore, if the conformal limit is reached for densities much larger than those considered here, the ultracompact region with radii 6–7 km can be populated [33]. In the opposite case of the early onset of the conformal limit (as discussed in Section 4), the radii will remain large, but small-mass regions can be populated if the stability criteria are modified by the slow conversion at the interface(s) between the phases. Another interesting new observation is that the change in the magnitude of the jump from 2SC to the CFL phase induces a special point on the M–R diagram at which all the curves meet in analogy to the case of single-phase transition; see Ref. [58]. The importance of studying this asymptotically large central density regime is phenomenologically relevant if the conversion between various quark and nuclear phases is slow compared to the characteristic timescale of oscillations, as discussed in Section 6. In this case, the stars on the descending branch of mass–central-density (and its counterpart on the M-R diagram) may be stable [17,18,19,20], contrary to the classical requirement $dM/d{\rho }_c>0$ for the branch to be stable, which in turn leads to higher multipole (beyond triplets) stars on the M–R diagram.