Phases and properties of color superconductors

Abstract

Cold and dense matter is expected to be in a color-superconducting state. Here we review two calculations, relevant for fundamental properties and applications of color superconductivity, respectively: the weak-coupling QCD calculation of the fermionic energy gap together with the magnetic screening masses of the gauge bosons, and the calculation of bulk viscosity from a non-leptonic electroweak process. These calculations are supplemented by a discussion of color superconductors with mismatched Fermi momenta, and they are embedded in the context of the state of the art by giving an overview of previous and ongoing work and future directions.

I Introduction

QCD matter at large baryon densities and small temperatures exhibits many phenomena known from ordinary, low-energy condensed matter physics. The most important one is Cooper pairing. Recall Cooper pairing in an ordinary superconductor: electrons immersed in a lattice of ions effectively attract each other, which leads to the formation of a Cooper pair condensate, which in turn induces an energy gap in the quasi-particle spectrum of the electrons. As a consequence, single electrons cannot be excited at low energies, which prevents them from undergoing scattering processes with the ion lattice while they are confined in a Cooper pair. Since electrons carry electric charge, this leads to the phenomenon of superconductivity. The underlying mechanism of Cooper pairing is very general. What we need is a sharp Fermi surface, i.e., a system of fermions in which the temperature is small compared to the Fermi energy. Then, if there is any attractive interaction between the fermions, even if it is arbitrarily small, Cooper pairing will occur. This is Cooper’s Theorem, and the underlying microscopic theory is BCS theory, after Bardeen, Cooper and Schrieffer PhysRev.108.1175. In an intuitive picture, Cooper pairing can be viewed as a way of fermions undergoing Bose-Einstein condensation. Since they are fermions, they cannot condense directly. So they pair, to form bosons, and then they condense. This picture has to be taken with some care, especially at weak coupling, where the size of a Cooper pair is much larger than the mean particle distance. However, experiments with ultra-cold atoms have shown that the transition to a strongly coupled system where fermion pairs are well-separated bosons is continuous. Therefore, the picture of Bose-Einstein condensation of Cooper pairs is, at least qualitatively, also applicable to weak coupling. Although best known in the context of electronic superconductivity, Cooper pairs do not necessarily carry electric charge. They are neutral for instance in ultra-cold fermionic gases or in ³He. In this case, the system is superfluid rather than superconducting. From a theoretical point of view, superconductors and superfluids are thus very similar. The difference lies in the charge of the Cooper pairs. Or, more formally speaking, superfluidity is associated with spontaneous symmetry breaking of a global symmetry (particle number conservation), while in a superconductor the broken group is local (the gauge group of electromagnetism or, as we shall see now, the color gauge group).

Due to the generality of Cooper’s Theorem, we may ask whether Cooper pairing also occurs in QCD. At high densities and low temperatures, there is a sharp Fermi sphere of quarks, which interact by the exchange of gluons. We know – for instance because of the existence of baryons – that there must be an attractive channel for this interaction. Therefore, we may apply Cooper’s Theorem and can expect the formation of Cooper pairs Barrois:1977xd; Frautschi:1978rz; Bailin:1983bm. Here the relevant interaction is the fundamental force of QCD. This is in contrast to electronic systems, where the attractive force is provided in a less direct way by the interaction with the ion lattice. Despite the common and very general mechanism, there are, of course, important differences between quark Cooper pairing and Cooper pairing of electrons or ³He atoms. For instance, in three-flavor quark matter, there are $N_{c}N_{f}=9$ fermion species, where $N_{c}$ is the number of colors and $N_{f}$ the number of flavors. As a consequence, there is a multitude of possible pairing patterns, i.e., if we want to describe a phase of dense quark matter we have to specify who pairs with who. Moreover, in principle all nine fermion species can have different Fermi momenta. Since weak-coupling Cooper pairing occurs in a small vicinity of a common Fermi surface, a splitting of Fermi momenta may break some of the Cooper pairs. Indeed, due to the strange quark mass and if we impose electric charge neutrality and equilibrium with respect to the weak interactions, (unpaired) quark matter at densities relevant for neutron stars does not have a single common Fermi surface. This complication can be neglected at ultra-high densities, where the quark chemical potential $\mu$ is much larger than all three quark masses, $\mu\gg m_{u}\simeq m_{d}\simeq m_{s}$, which can therefore be neglected. In this case, (unpaired) three-flavor quark matter has equal densities of up, down and strange quarks, $n_{u}=n_{d}=n_{s}$. It is thus electrically neutral and in electroweak equilibrium without the presence of any electrons or positrons, and Cooper pairing is undisrupted.

Quarks carry color charge, electric charge, and baryon number, and thus we may ask whether dense quark matter is a color superconductor, an electronic superconductor, and/or a superfluid. The answer depends on the pairing pattern and is best discussed in terms of the broken symmetries. At ultra-high densities, the color-flavor locked (CFL) phase Alford:1997zt; Alford:1998mk is energetically favored. The reason is that it is the only phase in which all nine quark species participate in pairing and thus it has the largest condensation energy. The CFL symmetry breaking pattern is

$$[SU(3)_{c}]\times\underbrace{SU(3)_{L}\times SU(3)_{R}}_{\displaystyle{\supset[U(1)_{Q}]}}\times U(1)_{B}\to\underbrace{SU(3)_{c+L+R}}_{\displaystyle{\supset[U(1)_{{\tilde{Q}}}]}}\times\mathbb{Z}_{2}\,.$$ (1)

Here, local symmetries are distinguished from global ones by square brackets. The left-hand side shows the symmetries of QCD: the color gauge group, the three-flavor chiral symmetry group (which is only approximate once quark masses are taken into account) and baryon number conservation, together with the QED gauge group, which is a subgroup of the flavor symmetry group. We have omitted the axial $U(1)_{A}$, which is broken on a quantum level at non-asymptotic densities. The right-hand side is the symmetry group under which the CFL state is invariant. This residual group locks transformations in color and flavor space. CFL is a color superconductor in the sense that the QCD gauge group is spontaneously broken, which manifests itself in nonzero Meissner masses of the gluons Rischke:2000ra; Schmitt:2003aa. CFL is a superfluid because baryon number conservation is broken. This manifests itself for instance in the formation of vortices if CFL is rotated (possibly in the interior of neutron stars) Forbes:2001gj; Balachandran:2005ev and in the existence of a Goldstone mode. This Goldstone mode is exactly massless because baryon number conservation is an exact symmetry of QCD. The electromagnetic gauge group $U(1)_{Q}$ survives in a “rotated” form, where the generator $\tilde{Q}$ is given by a linear combination of the original generator $Q$ and one of the generators of the color gauge group $SU(3)_{c}$. As a consequence, one gauge boson, a mixture of a gluon and the photon, remains massless.

The symmetry breaking pattern (1) also shows that the chiral symmetry group $SU(3)_{L}\times SU(3)_{R}$ is spontaneously broken. The Cooper pair condensate breaks it down to the same residual group as a chiral condensate does in low-density QCD. In the asymptotic limit $\mu\to\infty$, where the chiral group is an exact symmetry of QCD, the symmetry breaking pattern indicates the existence of eight Goldstone modes. At moderate densities, these Goldstone modes acquire a (small) mass since non-negligible quark masses break chiral symmetry explicitly. Due to the gap in the spectrum of all quasi-fermions, the pseudo-Goldstone modes from chiral symmetry breaking and the exact superfluid Goldstone mode dominate the transport properties of CFL, and one can write down an effective low-energy theory Son:1999cm; Bedaque:2001je, very similar to chiral perturbation theory in the QCD vacuum.

This review article is structured as follows. In Sec. II, we shall go through the weak-coupling calculation of the color-superconducting gap, the resulting condensation energy, and the Meissner masses of the gauge bosons. These calculations provide rigorous QCD results which, due to asymptotic freedom, apply to asymptotically large densities. We shall mostly restrict ourselves to CFL, but we will mention other possible phases and how they differ from CFL. In Sec. III we discuss – in a less detailed way – the fate of the CFL phase as we go down in density. We present another detailed and exemplary calculation in Sec. IV, where we calculate the bulk viscosity in the 2SC phase, making use of the quark propagators in the presence of Cooper pairing discussed in Sec. II. The bulk viscosity is one example of a transport coefficient, and Sec. IV points out the typical features encountered in transport of (partially) paired fermions. This section also contains a comparison to the bulk viscosity of CFL and a brief general discussion of transport in neutron stars, focusing on the effect of color superconductivity. By concentrating on the pedagogical presentation of two detailed calculations many aspects of color superconductivity can only be touched upon very briefly. More details can be found in the references and in the already existing reviews, which are recommended for further reading Rajagopal:2000wf; Buballa:2003qv; Shovkovy:2004me; Shovkovy:2005fy; Alford:2007xm; Wang:2009xf; Schmitt:2010pn; Fukushima:2010bq; Pang:2011qh; Anglani:2013gfu; Mannarelli:2014jsa.

II Cooper pairing in QCD

II.1 Quark propagator

A very useful formalism for many-body systems with spontaneous symmetry breaking is the 2-particle-irreducible (2PI) formalism, also called Cornwall-Jackiw-Tomboulis (CJT) formalism Cornwall:1974vz; Baym:1962sx; Luttinger:1960ua. The starting point is an effective action $\Gamma$, which is a functional of the quark and gluon propagators $S$ and $D$,

$$\Gamma[S,D]=\frac{1}{2}{\rm Tr}\ln S^{-1}-\frac{1}{2}{\rm Tr}(1-S_{0}^{-1}S)-\frac{1}{2}{\rm Tr}\ln D^{-1}+\frac{1}{2}{\rm Tr}(1-D_{0}^{-1}D)+\Gamma_{2}[S,D]\,,$$ (2)

where $S_{0}$ and $D_{0}$ are the free (non-interacting) propagators. The factor 1/2 in front of the fermionic terms is due to the use of Nambu-Gorkov propagators, which is necessary in the presence of Cooper pairing, see below. The contribution $\Gamma_{2}[S,D]$ contains all 2-particle-irreducible diagrams, i.e., diagrams that remain connected after cutting one or two lines. Truncating at two-loop order, the relevant diagrams are shown in Fig. 1.

Figure 1: Diagrams for the 2-loop truncation of the 2PI effective action. Solid (curly) lines correspond to the quark (gluon) propagator $S$ ($D$). The pressure is dominated by the left diagram because the gluon contributions are suppressed at small temperatures. Also for the other quantities calculated here the middle and right diagrams are not relevant.

The full propagators are obtained by the stationarity equations, which are obtained by taking the functional derivatives of the effective action,

$$0=\frac{\delta\Gamma}{\delta S}=\frac{\delta\Gamma}{\delta D}\,,$$ (3)

which yields

$$S^{-1} = S_{0}^{-1}+\Sigma\,,$$ (4a)

$$D^{-1} = D_{0}^{-1}+\Pi\,,$$ (4b)

where

$$\Sigma=2\frac{\delta\Gamma_{2}}{\delta S}\,,\qquad\Pi=-2\frac{\delta\Gamma_{2}}{\delta D}$$ (5)

are the quark and gluon self-energies, respectively, see Fig. 2.

Figure 2: The quark self-energy $\Sigma$ (left diagram) is obtained by cutting a fermion line in the left diagram of Fig. 1. The gluon self-energy $\Pi$ is obtained by cutting a gluon line in the diagrams of Fig. 1. We will only be interested in the contribution to the gluon self-energy from the quark loop (right diagram).

The quark self-energy $\Sigma$ will be relevant for the QCD gap equation, see Sec. II.2. The gluon self-energy $\Pi$ – more precisely, the quark loop contribution to it – will be needed in two instances: for the gap equation it is sufficient to use the so-called Hard-Dense-Loop (HDL) approximation for $\Pi$, where the quark loop contains ungapped quarks, while for the calculation of the screening masses from $\Pi$ the gapped quark propagator is needed, see Sec. II.4.

The grand-canonical potential (or free energy density) $\Omega$ is obtained by evaluating the effective action at the stationary point. With $\Gamma_{2}=\frac{1}{4}{\rm Tr}(\Sigma S)$ we obtain

$$\Omega\simeq-\frac{1}{2}{\rm Tr}\ln S^{-1}+\frac{1}{4}{\rm Tr}(1-S_{0}^{-1}S)\,.$$ (6)

Here we have neglected the contribution of the gluons because we are only interested in small temperatures $T\ll\mu$, where their contribution is negligible.

With the formalism at hand, we now include Cooper pairing. Cooper pairing gives rise to a nonzero expectation value of two fermion fields, so naively we expect a nonzero $\langle\psi\psi\rangle$, where $\psi$ is the quark spinor. However, the product $\psi\psi$ is not defined (think of $\psi$ as a column vector in Dirac space). This is in contrast to the product $\bar{\psi}\psi$, which appears in the chiral condensate $\langle\bar{\psi}\psi\rangle$ and which can be understood as the product of a row vector with a column vector, yielding a scalar. Instead, Cooper pairing is described by a condensate of the form $\langle\psi\bar{\psi}_{C}\rangle$, where $\psi_{C}=C\bar{\psi}^{T}$ with $C=i\gamma^{2}\gamma^{0}$ is the charge-conjugate spinor. (In $\bar{\psi}_{C}$, first charge-conjugate, then take the Hermitian conjugate and then multiply by $\gamma^{0}$.) It is therefore convenient to work in an extended spinor space, where fermions and charge-conjugate fermions are subsumed in the so-called Nambu-Gorkov spinor

$$\Psi=\left(\begin{array}[]{c}\psi\\ \psi_{C}\end{array}\right)\,.$$ (7)

As a consequence, the fermion propagators $S_{0}^{-1}$ and $S$ acquire an additional $2\times 2$ structure. Together with Dirac, color, and flavor structure, they are therefore $72\times 72$ matrices ($2\times 4N_{c}N_{f}=72$). The free inverse propagator in Nambu-Gorkov space is

$$S_{0}^{-1}=\left(\begin{array}[]{cc}[G_{0}^{+}]^{-1}&0\\ 0&[G_{0}^{-}]^{-1}\end{array}\right)\,,$$ (8)

where the free inverse fermion propagator and the free inverse charge-conjugate fermion propagator in momentum space are

$$[G_{0}^{\pm}]^{-1}=\gamma^{\mu}K_{\mu}\pm\mu\gamma^{0}=\sum_{e=\pm}[k_{0}\pm(\mu-ek)]\gamma^{0}\Lambda_{k}^{\pm e}\,,$$ (9)

with the four-momentum $K^{\mu}=(k_{0},\bm{k})$, and $k=|\bm{k}|$. The temporal component is given by the fermionic Matsubara frequencies, $k_{0}=-i\omega_{n}$, with $\omega_{n}=(2n+1)\pi T$, where $n\in\mathbb{Z}$. We have restricted ourselves for simplicity to massless quarks (having in mind asymptotically large densities), and we have introduced the energy projectors

$$\Lambda_{k}^{e}=\frac{1}{2}(1+e\gamma^{0}\bm{\gamma}\cdot\hat{\bm{k}})\,,$$ (10)

which are orthogonal, $\Lambda_{k}^{+}\Lambda_{k}^{-}=0$, and complete, $\Lambda_{k}^{+}+\Lambda_{k}^{-}=1$. For the following, it is useful to note that $\gamma^{0}\Lambda_{k}^{e}=\Lambda_{k}^{-e}\gamma^{0}$. Writing the propagator in terms of projectors is very useful for inversion. For instance, from Eq. (9) we immediately obtain

$$G_{0}^{\pm}=\sum_{e}\frac{\Lambda_{k}^{\pm e}\gamma^{0}}{k_{0}\pm(\mu-ek)}\,.$$ (11)

For now, by defining the inverse Nambu-Gorkov propagator (8) we have simply doubled the degrees of freedom since fermions and charge-conjugate fermions remain uncoupled. They couple in the self-energy, which we write as

$$\Sigma=\left(\begin{array}[]{cc}\Sigma^{+}&\Phi^{-}\\ \Phi^{+}&\Sigma^{-}\end{array}\right)\,.$$ (12)

The diagonal elements can be approximated by Gerhold:2005uu

$$\Sigma^{\pm}\simeq\gamma^{0}\Lambda_{k}^{\pm}\frac{g^{2}}{18\pi^{2}}k_{0}\ln\frac{48e^{2}m_{g}^{2}}{\pi^{2}k_{0}^{2}}\,,$$ (13)

where $e$ is Euler’s constant, $g$ is the strong coupling constant [$\alpha_{s}=g^{2}/(4\pi)$], and

$$m_{g}^{2}\equiv N_{f}\frac{g^{2}\mu^{2}}{6\pi^{2}}$$ (14)

is the effective gluon mass squared in a dense QCD plasma with $N_{f}$ quark flavors and $\mu\gg T$. The off-diagonal elements of the self-energy $\Sigma$ describe Cooper pairing. We will refer to $\Phi^{+}$ as the gap matrix; $\Phi^{-}$ is related to $\Phi^{+}$ via $\Phi^{-}=\gamma^{0}(\Phi^{+})^{\dagger}\gamma^{0}$. The gap matrix is a matrix in color, flavor, and Dirac space and determines the pairing pattern. We shall continue the calculation with the ansatz

$$\Phi^{+}=\Delta{\cal M}\gamma^{5}\,.$$ (15)

The gap function $\Delta$ depends on the four-momentum $K$ and will be determined by solving the gap equation. The Dirac structure $\gamma^{5}$ ensures that Cooper pairs are made of quarks with the same chirality and have total spin zero. In general, the Dirac structure can be more complicated, for instance if pairing occurs in the spin-1 channel. For a discussion of a more general Dirac structure see Refs. Bailin:1983bm; Pisarski:1999av. In our ansatz Dirac and color/flavor structures factorize, such that ${\cal M}$ is a $9\times 9$ matrix in color/flavor space. In order to discuss the structure of ${\cal M}$ we consider the irreducible representations of diquarks regarding color and flavor symmetries,

$$SU(3)_{c}:\qquad[3]_{c}\otimes[3]_{c} = [\bar{3}]_{c}^{a}\oplus[6]^{s}_{c}\,,$$ (16a)

$$SU(3)_{f}:\qquad[3]_{f}\otimes[3]_{f} = [\bar{3}]_{f}^{a}\oplus[6]^{s}_{f}\,.$$ (16b)

Here, the flavor sector ‘$f$’ represents both left-handed and right-handed sectors. The anti-symmetric anti-triplet channel in color space $[\bar{3}]_{c}^{a}$ is attractive. Since the overall wave function of a quark Cooper pair has to be anti-symmetric with respect to exchange of the two fermions, this determines the flavor channel: if pairing occurs in the anti-symmetric spin-0 channel, the flavor channel has to be anti-symmetric as well, and thus ${\cal M}\in[\bar{3}]_{c}^{a}\otimes[\bar{3}]_{f}^{a}$. The general form of the color/flavor structure can thus be written as

$${\cal M}=\phi_{AB}J_{A}\otimes I_{B}\,,\qquad{\cal M}^{ij}_{\alpha\beta}=-\phi_{AB}\epsilon_{A\alpha\beta}\epsilon_{Bij}\,,$$ (17)

where $\{J_{1},J_{2},J_{3}\}$ with $(J_{A})_{\alpha\beta}=-i\epsilon_{A\alpha\beta}$ and $\{I_{1},I_{2},I_{3}\}$ with $(I_{B})_{ij}=-i\epsilon_{Bij}$ are bases of $[\bar{3}]_{c}^{a}$ and $[\bar{3}]_{f}^{a}$, respectively. The $3\times 3$ matrix $\phi$ determines the pairing pattern and thus the particular color-superconducting phase. We shall mostly be concerned with CFL, but also mention the so-called 2SC phase,

$${\rm CFL:}\quad\phi_{AB}=\delta_{AB}\,,\qquad{\rm 2SC:}\quad\phi_{AB}=\delta_{A3}\delta_{B3}\,.$$ (18)

In the 2SC phase, quarks of one color, say blue, remain unpaired, as well as the quarks of one flavor. If all quark masses can be neglected (and if we consider pure QCD, i.e., if the electric charges of the quarks play no role) there is no difference between the flavors and we can randomly pick one quark flavor that remains unpaired. However, the 2SC phase becomes particularly interesting if the strange quark mass is taken into account. In that case, the 2SC phase usually refers to the phase where strange quarks remain unpaired. As we shall see below, at asymptotically large densities the CFL phase is favored. (In the presence of a large magnetic field the 2SC phase may become favored even in the massless limit Haber:2017oqb.) One can check that the CFL order parameter in Eq. (18) gives rise to the symmetry breaking pattern (1).

We may now compute the full quark propagator. We shall do so for the CFL phase and neglect the diagonal elements of the self-energy, $\Sigma^{\pm}\simeq 0$. Their effect is not very crucial for the understanding of the main physics, and neglecting them renders the calculation simpler. At the end, for completeness we will reinstate their (subleading) effect on the superconducting gap. With this simplification, putting together equations (4a), (8), and (12), we obtain the inverse propagator

$$S^{-1}=\left(\begin{array}[]{cc}[G_{0}^{+}]^{-1}&\Phi^{-}\\ \Phi^{+}&[G_{0}^{-}]^{-1}\end{array}\right)\,.$$ (19)

Inversion of this matrix yields

$$S=\left(\begin{array}[]{cc}G^{+}&F^{-}\\ F^{+}&G^{-}\end{array}\right)\,,$$ (20)

with

$$G^{\pm} = \left([G_{0}^{\pm}]^{-1}-\Phi^{\mp}G_{0}^{\mp}\Phi^{\pm}\right)^{-1}\,,$$ (21a)

$$F^{\pm} = -G_{0}^{\mp}\Phi^{\pm}G^{\pm}\,.$$ (21b)

The off-diagonal elements in Nambu-Gorkov space $F^{\pm}$ are called anomalous propagators. They are only nonzero in the presence of Cooper pairing and describe (upper sign) the propagation of a charge-conjugate fermion which turns into a propagating fermion through the Cooper pair condensate. This is a manifestation of the spontaneous breaking of baryon number conservation: fermions can be created and annihilated by the condensate. In order to perform the inversion in Eq. (21a) it is useful to write the matrix ${\cal M}^{\dagger}{\cal M}$ (for CFL, ${\cal M}^{\dagger}={\cal M}$) in its spectral representation,

$${\cal M}^{2}=\lambda_{1}{\cal P}_{1}+\lambda_{2}{\cal P}_{2}\,,$$ (22)

where $\lambda_{1}=1$ and $\lambda_{2}=4$ are the eigenvalues of ${\cal M}^{2}$, and ${\cal P}_{1}$ and ${\cal P}_{2}$ are the projectors onto the corresponding eigenspaces,

$${\cal P}_{1}=-\frac{{\cal M}^{2}-4}{3}\,,\qquad{\cal P}_{2}=\frac{{\cal M}^{2}-1}{3}\,.$$ (23)

The degeneracies of the eigenvalues are ${\rm Tr}({\cal P}_{1})=8$ and ${\rm Tr}({\cal P}_{2})=1$. In terms of color and flavor indices ($\alpha,\beta$ and $i,j$, respectively) we have

$${\cal M}_{\alpha\beta}^{ij}=\delta_{\alpha}^{j}\delta_{\beta}^{i}-\delta_{\alpha}^{i}\delta_{\beta}^{j}\quad\Rightarrow\qquad({\cal M}^{2})_{\alpha\beta}^{ij}=\delta_{\alpha}^{i}\delta_{\beta}^{j}+\delta_{\alpha\beta}\delta^{ij}\,,$$ (24)

which yields

$${\cal P}_{1}=\delta_{\alpha\beta}\delta^{ij}-\frac{1}{3}\delta_{\alpha}^{i}\delta_{\beta}^{j}\,,\qquad{\cal P}_{2}=\frac{1}{3}\delta_{\alpha}^{i}\delta_{\beta}^{j}\,.$$ (25)

From these expressions one finds ${\cal M}^{4}=5{\cal M}^{2}-4$, which can be used to show that the projectors are orthogonal, ${\cal P}_{1}{\cal P}_{2}=0$. They are, obviously, also complete, ${\cal P}_{1}+{\cal P}_{2}=1$. With the help of these projectors we compute from Eq. (21)

$$G^{\pm} = \sum_{e=\pm}\sum_{r=1,2}\frac{[k_{0}\mp(\mu-ek)]{\cal P}_{r}\gamma^{0}\Lambda_{k}^{\mp e}}{k_{0}^{2}-(\epsilon_{k,r}^{e})^{2}}\,,$$ (26a)

$$F^{\pm} = \sum_{e=\pm}\sum_{r=1,2}\frac{\Phi^{\pm}{\cal P}_{r}\Lambda_{k}^{\mp e}}{k_{0}^{2}-(\epsilon_{k,r}^{e})^{2}}\,,$$ (26b)

where

$$\epsilon_{k,r}^{e}=\sqrt{(\mu-ek)^{2}+\lambda_{r}\Delta^{2}}\,.$$ (27)

This result shows that the propagator has poles at $k_{0}=\pm\epsilon_{k,r}^{e}$. These are the quasi-particle and quasi-antiparticle excitations (upper sign) and quasi-hole and quasi-antihole excitations (lower sign), see Fig. 3.

Figure 3: Dispersion relations with (solid) and without (dashed) Cooper pairing. In the absence of Cooper pairing, quasi-fermions can be excited with infinitesimally small energy at the Fermi surface. If Cooper pairs have formed, an energy $2\Delta$ is necessary to excite quasi-fermions. These quasi-fermions are (momentum-dependent) mixtures of particles and holes. Also anti-particle (uppermost branch) and anti-hole (lowermost branch) excitations are shown for completeness.

The structure of these excitations is solely determined by the pairing pattern, i.e., the spectrum of the matrix ${\cal M}^{2}$: in the CFL phase, 8 quasi-particles are gapped with gap $\Delta$, while one quasi-particle is gapped with gap $2\Delta$. From the matrix ${\cal M}$ one can also read off the pairing pattern: if we label the color directions by $r$, $g$, $b$ and the flavor directions with $u$, $d$, $s$, then we have pairing of $rd$ with $gu$, of $bu$ with $rs$, of $bd$ with $gs$, and pairing among the three quark species $ru$, $gd$, $bs$. The calculation for the 2SC phase is completely analogous and yields 4 gapped quasi-particles and 5 ungapped quasi-particles. Here, $ru$ pairs with $gd$, $rd$ with $gu$, while $bu$, $bd$, $bs$, $rs$, $gs$ remain unpaired.

II.2 Solving the gap equation

To compute the gap function $\Delta(K)$ we go back to the quark-self energy $\Sigma$, which is given by the left diagram in Fig. 2. This diagram contains the full propagator in Nambu-Gorkov space, i.e., it has to be understood as a $2\times 2$ matrix, whose four elements are given by the quark propagators $G^{\pm}$ and $F^{\pm}$, see Eq. (20). On the other hand, the self-energy is given by the gap matrix according to Eq. (12). Therefore, the lower off-diagonal element in Nambu-Gorkov space relates the gap matrix to the self-energy diagram. This is the gap equation, which is shown diagrammatically in Fig. 4 (equivalently, we could have considered the upper off-diagonal element).

Figure 4: Diagrammatic representation of the QCD gap equation. The right-hand side is the (lower) off-diagonal element in Nambu-Gorkov space of the left diagram in Fig. 2, containing the anomalous propagator $F^{+}=-G_{0}^{-}\Phi^{+}G^{+}$. The vertices are given by $g\gamma^{\mu}T_{a}^{T}$ and $g\gamma^{\nu}T_{b}$, whose indices are contracted by the gluon propagator $D_{\mu\nu}^{ab}$.

The algebraic version of the gap equation is

$$\Phi^{+}(K)=g^{2}\frac{T}{V}\sum_{Q}\gamma^{\mu}T_{a}^{T}F^{+}(Q)\gamma^{\nu}T_{b}D_{\mu\nu}^{ab}(P)\,,$$ (28)

where $V$ is the volume of the system, where $T_{a}=\lambda_{a}/2$ with the Gell-Mann matrices $\lambda_{a}$ ($a=1,\ldots,8$), and where we have denoted the gluon four-momentum by $P\equiv K-Q$. For our purposes we may assume the gluon propagator to be diagonal in color space, $D_{\mu\nu}^{ab}=\delta^{ab}D_{\mu\nu}$, and employ the HDL approximation. In this approximation, the gluon-self energy is computed from a quark loop where the quark momentum is ‘hard’, i.e., of the order of the chemical potential, and much larger than the ‘soft’ gluon energy and momenta. In particular, in this approximation, the effect of Cooper pairing is not taken into account. In other words, in computing the color-superconducting gap, the backreaction of the gap on the gluon propagator, which in turn affects the gap, can be neglected up to the order we will be working Rischke:2001py.

In this approximation the gluon self-energy $\Pi$ is (4-)transverse to the gluon momentum $P$, as in abelian gauge theories. This allows us to decompose the self-energy into (3-)longitudinal and (3-)transverse components with respect to $\bm{q}$,

$$\Pi_{\mu\nu}={\cal F}P_{L,\mu\nu}+{\cal G}P_{T,\mu\nu}\,,$$ (29)

where transverse and longitudinal projectors are defined as

$$P_{T,00} = P_{T,0i}=P_{T,i0}=0\,,\qquad P_{T,ij}=\delta_{ij}-\hat{p}_{i}\hat{p}_{j}\,,$$ (30a)

$$P_{L,\mu\nu} = \frac{P_{\mu}P_{\nu}}{P^{2}}-g_{\mu\nu}-P_{T,\mu\nu}\,.$$ (30b)

The components ${\cal F}$ and ${\cal G}$ can be computed from the self-energy by appropriate projections,

$${\cal F}=\frac{P^{2}}{p^{2}}\Pi_{00}\,,\qquad{\cal G}=\frac{1}{2}(\delta_{ij}-\hat{p}_{i}\hat{p}_{j})\Pi_{ji}\,.$$ (31)

In the HDL approximation Bellac:2011kqa,

$${\cal F}(P) = -3m_{g}^{2}\frac{P^{2}}{p^{2}}\left(1-\frac{p_{0}}{2p}\ln\frac{p_{0}+p}{p_{0}-p}\right)\,,$$ (32a)

$${\cal G}(P) = \frac{3m_{g}^{2}p_{0}}{2p}\left(\frac{p_{0}}{p}-\frac{P^{2}}{2p^{2}}\ln\frac{p_{0}+p}{p_{0}-p}\right)\,,$$ (32b)

with the effective gluon mass $m_{g}$ from Eq. (14). In Coulomb gauge, the gluon propagator in the HDL approximation is

$$D_{00}=D_{L}\,\qquad D_{0i}=0\,,\qquad D_{ij}=(\delta_{ij}-\hat{p}_{i}\hat{p}_{j})D_{T}\,,$$ (33)

where the longitudinal and transverse components are

$$D_{L}(P)=\frac{P^{2}}{p^{2}}\frac{1}{{\cal F}(P)-P^{2}}\,,\qquad D_{T}(P)=\frac{1}{{\cal G}(P)-P^{2}}\,.$$ (34)

The physical meaning of this form of the propagator is best understood by computing the corresponding spectral densities, which are given by the imaginary part of the retarded propagators,

$$\rho_{L,T}=\frac{1}{\pi}\lim_{\epsilon\to 0}{\rm Im}\,D_{L,T}(p_{0}+i\epsilon,\bm{p})\,,$$ (35)

where $\epsilon>0$. With Eqs. (32) and (34) one computes

$$\rho_{L} = \frac{\omega_{L}(\omega_{L}^{2}-p^{2})\delta(p_{0}-\omega_{L})}{p^{2}(3m_{g}^{2}+p^{2}-\omega_{L}^{2})}-\frac{1}{\pi}\frac{P^{2}}{p^{2}}\frac{\Theta(p-p_{0}){\rm Im}\,{\cal F}}{({\rm Re}\,{\cal F}-P^{2})^{2}+({\rm Im}\,{\cal F})^{2}}\,,$$ (36a)

$$\rho_{T} = \frac{\omega_{T}(\omega_{T}^{2}-p^{2})\delta(p_{0}-\omega_{T})}{3m_{g}^{2}\omega_{T}^{2}-(\omega_{T}^{2}-q^{2})^{2}}-\frac{1}{\pi}\frac{\Theta(p-p_{0}){\rm Im}\,{\cal G}}{({\rm Re}\,{\cal G}-P^{2})^{2}+({\rm Im}\,{\cal G})^{2}}\,,$$ (36b)

where $\omega_{L}$ and $\omega_{T}$ are the poles of the propagator, i.e., the numerical solutions for $p_{0}$ of $\Pi_{00}-p^{2}=0$ and ${\cal G}-P^{2}=0$, respectively. We see that the gluon propagator exhibits two quasi-particle excitations in the timelike regime, $P^{2}>0$, given by the delta peak in the spectral density. These are the longitudinal and transverse plasmon modes. In addition, the spectral density is nonzero in the entire spacelike regime, $P^{2}<0$. These are Landau-damped gluons, unstable excitations due to the scattering with the quarks in the plasma, which is a well-known phenomenon also for photons in a QED plasma. In fact, in the given HDL approximation, the diagonal elements of the gluon propagator are identical to the photon propagator in an electromagnetic plasma.

Coming back to the gap equation (28), we insert the anomalous propagator (26b) and neglect the antiparticle contribution, i.e., we drop the term with $e=-$ and from now on denote $\epsilon_{k,r}\equiv\epsilon_{k,r}^{+}$. Moreover, we use the form of the gluon propagator from Eq. (33), multiply both sides of the gap equation by $\gamma^{5}{\cal M}\Lambda_{k}^{+}$ from the right and take the trace over color, flavor, and Dirac space. With the color/flavor traces

$${\rm Tr}[T_{a}^{T}{\cal M}{\cal P}_{1}T_{a}{\cal M}]=2{\rm Tr}[T_{a}^{T}{\cal M}{\cal P}_{2}T_{a}{\cal M}]=-\frac{16}{3}\,,$$ (37)

and the Dirac traces

$${\rm Tr}[\gamma^{0}\gamma^{5}\Lambda_{q}^{-}\gamma^{0}\gamma^{5}\Lambda_{k}^{+}] = -(1+\hat{\bm{q}}\cdot\hat{\bm{k}})\,,$$ (38a)

$${\rm Tr}[\gamma^{i}\gamma^{5}\Lambda_{q}^{-}\gamma^{j}\gamma^{5}\Lambda_{k}^{+}] = \delta_{ij}(1-\hat{\bm{q}}\cdot\hat{\bm{k}})+\hat{q}_{i}\hat{k}_{j}+\hat{q}_{j}\hat{k}_{i}\,,$$ (38b)

we find

$$\Delta(K) = \frac{g^{2}}{3}\frac{T}{V}\sum_{Q}\left(\frac{2}{3}\frac{\Delta(Q)}{q_{0}^{2}-\epsilon_{q,1}^{2}}+\frac{1}{3}\frac{\Delta(Q)}{q_{0}^{2}-\epsilon_{q,2}^{2}}\right)\left[(1+\hat{\bm{q}}\cdot\hat{\bm{k}})D_{L}(P)-2(1-\hat{\bm{p}}\cdot\hat{\bm{q}}\,\hat{\bm{p}}\cdot\hat{\bm{k}})D_{T}(P)\right]\,.$$ (39)

The sum over 4-momenta $Q$ contains the sum over Matsubara frequencies and the sum over 3-momenta, which turns into an integral in the thermodynamic limit,

$$\frac{T}{V}\sum_{Q}\to T\sum_{q_{0}}\int\frac{d^{3}\bm{q}}{(2\pi)^{3}}\,.$$ (40)

It is instructive to discuss the gap equation (39) first in the simple case of a pointlike interaction between fermions. So, let us (unrealistically) assume that quarks interact by the exchange of a heavy scalar field with mass $M$ such that we can approximate $D_{L}\simeq-1/M^{2}$ and $D_{T}=0$. In this case, the sum over the fermionic Matsubara frequencies is easily performed with the help of

$$T\sum_{q_{0}}\frac{\Delta(Q)}{q_{0}^{2}-\epsilon_{q}^{2}}=-\frac{\Delta_{q}}{2\epsilon_{q}}\tanh\frac{\epsilon_{q}}{2T}\,,$$ (41)

where $\Delta_{q}\equiv\Delta(\epsilon_{q},\bm{q})$. Moreover, we can simply ignore the energy and momentum dependence of the gap function and thus divide the whole equation by $\Delta$. ($\Delta=0$ is obviously a solution to the gap equation for all temperatures. Eventually one has to check that the phase with a nonzero $\Delta$ has lower free energy than the unpaired phase.) We perform the (trivial) angular integral and restrict the integral over momenta $q$ to a small vicinity of the Fermi surface, $q\in[\mu-\delta,\mu+\delta]$, where $\delta$ is much larger than the gap but much smaller than the Fermi momentum, $\Delta\ll\delta\ll\mu$. This allows us to approximate the integration measure by $dq\,q^{2}\simeq\mu^{2}dq$. Finally, with the new integration variable $\xi=q-\mu$ we obtain

$$1\simeq\frac{g^{2}}{c}\int_{0}^{\delta}d\xi\left(\frac{2}{3\epsilon_{q,1}}\tanh\frac{\epsilon_{q,1}}{2T}+\frac{1}{3\epsilon_{q,2}}\tanh\frac{\epsilon_{q,2}}{2T}\right)\,,$$ (42)

with the dimensionless constant $c\equiv 6M^{2}\pi^{2}/\mu^{2}$. At zero temperature, we have $\tanh\frac{\epsilon_{q,r}}{2T}=1$ and we can perform the resulting integral analytically. Denoting the zero-temperature gap by $\Delta_{0}$ and using $\delta\gg\Delta_{0}$, we find

$$\mbox{pointlike interaction:}\qquad\Delta_{0}=2\delta\cdot 2^{-1/3}\exp\left(-\frac{c}{g^{2}}\right)\,.$$ (43)

This detour has shown us, firstly, how to solve the gap equation in the BCS scenario, where there are no long-range interactions. Notice that the effective dimensional reduction of the momentum integral was crucial and can thus be considered as the formal reason for the instability towards Cooper pair condensation. Secondly, the result shows the BCS behavior of the gap: at weak coupling – where this calculation is valid – the gap is exponentially suppressed. The result is non-perturbative in the sense that there is no Taylor expansion of $\Delta_{0}$ for small $g$. By solving the gap equation, we have implicitly resummed infinitely many diagrams since the quark self energy, from which the gap is computed, contains the full quark propagator, which in turn contains the gap etc. Thirdly, we see that the two-gap structure of the CFL phase creates a factor $2^{-1/3}$ in the gap. Below, we shall quote the general form of this factor, assuming a two-gap structure but keeping the degeneracies and prefactors of the gap in the dispersion relation general.

One can use Eq. (42) also to compute the temperature dependence of the gap. This has to be done numerically, and one finds a monotonically decreasing gap which goes to zero at a certain temperature, the critical temperature $T_{c}$. This temperature can be calculated analytically. To this end, we set $\Delta=0$ in the excitation energies, such that $\epsilon_{q,1}=\epsilon_{q,2}=\xi$. As a consequence, the two-gap structure disappears and we have

$$1=\frac{g^{2}}{c}\int_{0}^{\delta}\frac{d\xi}{\xi}\tanh\frac{\xi}{2T_{c}}\,.$$ (44)

With the new integration variable $z=\xi/(2T_{c})$ we compute via partial integration

$$\frac{c}{g^{2}}=\ln z\,\tanh z\Big|_{0}^{\delta/(2T_{c})}-\int_{0}^{\delta/(2T_{c})}dz\frac{\ln z}{\cosh^{2}z}\simeq\ln\frac{\delta}{2T_{c}}+\gamma-\ln\frac{\pi}{4}\,,$$ (45)

where $\gamma\simeq 0.577$ is the Euler-Mascheroni constant, and where we have used $\delta\gg T_{c}$, such that the upper boundary of the integral can be approximated by infinity. This is a consistent assumption since, as we shall see now, the critical temperature is of the same order as the zero-temperature gap. Solving Eq. (45) for $T_{c}$ yields

$$T_{c}=\frac{e^{\gamma}}{\pi}2\delta e^{-c/g^{2}}=\frac{e^{\gamma}}{\pi}2^{1/3}\Delta_{0}\simeq 2^{1/3}\cdot 0.567\Delta_{0}\,.$$ (46)

We see that the nontrivial factor due to the two-gap structure, in this case $2^{1/3}$, appears in the relation between $T_{c}$ and $\Delta_{0}$, but not if $T_{c}$ is expressed directly in terms of the coupling constant.

Let us return to QCD and realistic gluons. For brevity we ignore the two-gap structure and set $\epsilon_{1,q}=\epsilon_{2,q}\equiv\epsilon_{q}$. Its effect on the QCD gap is the same as just shown for the BCS case. We may thus simply reinstate the resulting factor in the final result. It is convenient to work with the spectral representation of the gluon propagator, based on the spectral densities (36). Details of this procedure are worked out in Ref. Pisarski:1999tv. It turns out that the gap equation is sensitive to gluon momenta with $p_{0}\ll p$. This allows us to simplify the gluon propagator from the beginning,

$$D_{L}(p) \simeq -\frac{1}{p^{2}+3m_{g}^{2}}\,,$$ (47a)

$$D_{T}(p_{0},p) \simeq \frac{\Theta(p-M_{g})}{p^{2}}+\frac{\Theta(M_{g}-p)p^{4}}{p^{6}+M_{g}^{4}p_{0}^{2}}\,,$$ (47b)

where $M_{g}^{2}\equiv 3\pi m_{g}^{2}/4$. This yields, after performing the Matsubara sum and the (now nontrivial) angular integral,

$$\Delta_{k}\simeq\frac{g^{2}}{24\pi^{2}}\int_{\mu-\delta}^{\mu+\delta}dq\frac{\Delta_{q}}{\epsilon_{q}}\left(\ln\frac{4\mu^{2}}{3m_{g}^{2}}+\ln\frac{4\mu^{2}}{M_{g}^{2}}+\frac{1}{3}\ln\frac{M_{g}^{2}}{|\epsilon_{q}^{2}-\epsilon_{k}^{2}|}\right)\tanh\frac{\epsilon_{q}}{2T}\,.$$ (48)

The three terms in parentheses arise, respectively, from static electric gluons (47a), non-static magnetic gluons [first term of (47b)] and almost static, Landau-damped magnetic gluons [second term in (47b)]. Combining the three logarithms, we obtain

$$\Delta_{k}\simeq\bar{g}^{2}\int_{0}^{\delta}d(q-\mu)\frac{\Delta_{q}}{\epsilon_{q}}\frac{1}{2}\ln\frac{b^{2}\mu^{2}}{|\epsilon_{q}^{2}-\epsilon_{k}^{2}|}\tanh\frac{\epsilon_{q}}{2T}\,,$$ (49)

where

$$\bar{g}\equiv\frac{g}{3\sqrt{2}\pi}\,,\qquad b\equiv 256\pi^{4}\left(\frac{2}{N_{f}g^{2}}\right)^{5/2}\,.$$ (50)

We can now approximate the logarithm by

$$\frac{1}{2}\ln\frac{b^{2}\mu^{2}}{|\epsilon_{q}^{2}-\epsilon_{k}^{2}|}\simeq\Theta(k-q)\ln\frac{b\mu}{\epsilon_{k}}+\Theta(q-k)\ln\frac{b\mu}{\epsilon_{q}}\,,$$ (51)

and introduce the logarithmic integration variable

$$y\equiv\bar{g}\ln\frac{2b\mu}{q-\mu+\epsilon_{q}}\,.$$ (52)

With

$$x\equiv\bar{g}\ln\frac{2b\mu}{k-\mu+\epsilon_{k}}\,,\qquad x^{*}\equiv\bar{g}\ln\frac{2b\mu}{\Delta_{0}}\,,\qquad x_{0}\equiv\bar{g}\ln\frac{b\mu}{\delta}\,,$$ (53)

where $\Delta_{0}$ is the zero-temperature value of the gap at the Fermi surface, the zero-temperature gap equation can be written as

$$\Delta(x)\simeq x\int_{x}^{x^{*}}dy\,\Delta(y)+\int_{x_{0}}^{x}dy\,y\,\Delta(y)\,.$$ (54)

Differentiating twice with respect to $x$, this becomes a differential equation, $\Delta^{\prime\prime}(x)=-\Delta(x)$, with boundary conditions $\Delta(x^{*})=\Delta_{0}$ and $\Delta^{\prime}(x^{*})=0$ (the gap peaks at the Fermi surface). This yields $\Delta(x)=\Delta_{0}\cos(x^{*}-x)$ and by reinserting this solution into the integral version of the gap equation (54) one finds $\cos x^{*}\simeq 0$ and thus $\Delta(x)\simeq\Delta_{0}\sin x$. The gap at the Fermi surface $\Delta_{0}$ is then found from $x^{*}\simeq 1$, which yields the final result

$$\Delta_{0}=2b\mu e^{b_{0}^{\prime}}e^{-\zeta}e^{-d}\exp\left(-\frac{3\pi^{2}}{\sqrt{2}g}\right)\,.$$ (55)

Here we have included various effects that we have omitted in our derivation. Our derivation would yield the same result with $b_{0}^{\prime}=d=\zeta=0$. The most important feature of the weak-coupling QCD gap is the leading behavior in $g$, which is different compared to a four-fermion interaction: comparing with Eq. (43) we see that the QCD gap is parametrically enhanced due to the different power of $g$ in the exponential Barrois:1979pv; Son:1998uk. This behavior arises from the exchange of almost static magnetic gluons [third term in Eq. (28)]. Due to the (one-loop) running of the coupling, $1/g^{2}$ increases logarithmically with $\mu$ at asymptotically large $\mu$. Therefore, $\exp(-{\rm const}/g)$ decreases more slowly than $1/\mu$ and as a consequence $\Delta_{0}$ increases with $\mu$ at asymptotically large $\mu$, while $\Delta_{0}/\mu$ decreases.

The origin and values of the prefactors of the gap (55) are as follows:

• The factor

  $$e^{b_{0}^{\prime}}=\exp\left(-\frac{\pi^{2}+4}{8}\right)\simeq 0.177$$ (56)

  arises if the quark self-energy (13) is carried through the calculation of the gap Wang:2001aq.

• From Eq. (43) we know that the two-gap structure of CFL generates a prefactor $2^{-1/3}$. This factor can be written more generally as $e^{-\zeta}$ with

  $$\zeta=\frac{1}{2}\frac{\langle{\rm Tr}({\cal P}_{1})\lambda_{1}\ln\lambda_{1}+{\rm Tr}({\cal P}_{2})\lambda_{2}\ln\lambda_{2}\rangle}{\langle{\rm Tr}({\cal P}_{1})\lambda_{1}+{\rm Tr}({\cal P}_{2})\lambda_{2}\rangle}\,.$$ (57)

  Using the same notation as above, $\lambda_{1}$ and $\lambda_{2}$ are the eigenvalues of ${\cal M}{\cal M}^{\dagger}$, while ${\cal P}_{1}$ and ${\cal P}_{2}$ are the projectors onto the respective eigenspaces, such that ${\rm Tr}({\cal P}_{1})$ and ${\rm Tr}({\cal P}_{2})$ are the degeneracies. For CFL, we verify that with $\lambda_{1}=1$, $\lambda_{2}=4$, ${\rm Tr}({\cal P}_{1})=8$ and ${\rm Tr}({\cal P}_{2})=1$ we obtain $e^{-\zeta}=2^{-1/3}$. If there is only a single isotropic gap, as in the 2SC phase, we have $\lambda_{1}=1$ and $\lambda_{2}=0$, and thus the prefactor becomes trivial, $e^{-\zeta}=1$. The general formula also allows for anisotropic gaps, with the angular brackets denoting angular average. This becomes relevant for spin-1 Cooper pairing, where ${\cal M}$ includes the spin structure of the gap matrix, such that the eigenvalues $\lambda_{r}$ contain the angular dependence of the gap in momentum space. In this case, $e^{-\zeta}$ becomes nontrivial even for a single gap Schmitt:2004et.

• Finally, the factor $e^{-d}$ becomes nontrivial only for spin-1 Cooper pairing. This factor, in contrast to the factor $e^{-\zeta}$, does depend on the details of the interaction. Depending on the particular spin-1 phase, one finds values of $d$ varying between $d=6$ for pairing of quarks with the same chirality and $d=4.5$ for opposite-chirality pairing Schafer:2000tw; Schmitt:2004et. This suppresses the weak-coupling gap in the spin-1 channel by 2-3 orders of magnitude compared to the spin-0 channel.

The zero-temperature gap (55) is a result from first principles, systematically collecting all contributions to $\ln(\Delta_{0}/\mu)$ of order $g^{-1}$, $\ln g$ , and $g^{0}$. (Higher order contributions arise for instance from terms of the HDL gluon propagator that we have neglected or terms beyond the HDL approximation, but also from going beyond the one-loop truncation of the quark self-energy that we have used to derive the gap equation, for instance from vertex corrections.) The result is valid for sufficiently large chemical potentials; it was estimated that it starts to lose its validity at around $\mu\lesssim 10^{5}\,{\rm GeV}$ Rajagopal:2000rs. This corresponds to chemical potentials many orders of magnitude larger than in the center of neutron stars. To get an idea of the size of the gap at more moderate densities we can extrapolate the QCD result down in density. According to the two-loop running of the coupling constant, we estimate that for neutron star interiors with $\mu\sim 400\,{\rm MeV}$ we have $\alpha_{s}\sim 1$ and thus $g\sim 3.5$. This yields a color-superconducting gap of $\Delta_{0}\sim 10\,{\rm MeV}$. Given that this is a bold extrapolation over many orders of magnitude, it is reassuring that phenomenological models such as the Nambu-Jona-Lasinio (NJL) model Nambu:1961tp; Nambu:1961fr give similar – although somewhat larger – results, $\Delta_{0}\sim(10-100)\,{\rm MeV}$ Buballa:2001gj; Steiner:2002gx; Abuki:2004zk; Blaschke:2005uj; Ruster:2005jc; Warringa:2006dk; Gholami:2024diy. These numbers suggest that the energy gap from quark Cooper pairing is significantly larger than the corresponding quantity from Cooper pairing of nucleons, which is calculated to be of the order of $1\,{\rm MeV}$ at most Wambach:1992ik; Schulze:1996zz; Fabrocini:2005lvb; Cao:2006gq.

It is of course desirable to understand the color-superconducting gap better at moderate densities, beyond comparing the weak-coupling extrapolation with phenomenological models. Possible approaches to capture strong-coupling effects are based on Dyson-Schwinger equations Nickel:2006vf; Marhauser:2006hy; Muller:2013pya; Muller:2016fdr; Murgana:2025dfa and the functional renormalization group Braun:2021uua; Gholami:2025afm. It has also been attempted to treat color superconductivity (or at least a simplified variant thereof) within the gauge-gravity correspondence, which makes the strong-coupling limit accessible; currently, however, only in theories (very) different from QCD Chen:2009kx; Basu:2011yg; BitaghsirFadafan:2018iqr; Faedo:2018fjw; Preau:2025ubr; CruzRojas:2025fzs. It was also pointed out that it is possible to obtain an upper limit for $\Delta_{0}$ from a completely different approach, independent of any model description and perturbative QCD calculation: by asking for the maximal change to the equation of state due to quark Cooper pairing that is thermodynamically consistent with the known low-density nuclear matter equation of state and with known astrophysical constraints from neutron star data, it was estimated within a Bayesian analysis that $\Delta_{0}\lesssim 450\,{\rm MeV}$ at $\mu=2.6\,{\rm GeV}$ Kurkela:2024xfh.

One can also use the gap equation (49) to compute the critical temperature for color superconductivity. Details of this calculation can be found for instance in Ref. Schmitt:2002sc. Expressed in terms of the zero-temperature gap, the result is the same as derived above for the pointlike interaction, see right-hand side of Eq. (46). Therefore, the critical temperature of CFL is $T_{c}=e^{\gamma}/\pi\,2^{1/3}\Delta_{0}\simeq 0.72\,\Delta_{0}$. With the above estimate of $\Delta_{0}\sim 10\,{\rm MeV}$ or larger, this suggests that the temperature in the interior of neutron stars is – except for the very early stages of the life of the star Prakash:2000jr – much smaller than the critical temperature of color superconductors. Even in neutron star mergers, where the temperature can be a few tens of MeV Hammond:2021vtv, it is conceivable that quark Cooper pairing survives.

At asymptotically large densities, color superconductors are type-I superconductors since the coherence length $\xi\sim 1/\Delta$ is much larger than the penetration depth $\lambda\simeq 1/(g\mu)$ and thus the Ginzburg-Landau parameter $\kappa=\lambda/\xi\ll 1/\sqrt{2}$. This is relevant for the critical temperature because it is known from electronic superconductors that gauge field fluctuations become important in the type-I regime. As a consequence, the finite-temperature phase transition becomes first order and the critical temperature receives a correction of order $g$. The corrected critical temperature $T_{c}^{*}$ for color superconductivity at weak coupling is Giannakis:2004xt,

$$T_{c}^{*}=T_{c}\left(1+\frac{\pi^{2}}{12\sqrt{2}}g\right)\,.$$ (58)

II.3 Condensation energy

We have seen that the gap equation has two solutions: the trivial one, corresponding to the unpaired phase, and the nontrivial one, corresponding to the formation of a Cooper pair condensate. To establish that the color-superconducting phase is indeed the ground state of QCD at large densities, we have to compute its free energy and check whether it is lower than the free energy of the unpaired phase. We shall go through this calculation in this section. In principle, there could also be other phases whose free energy competes with the unpaired and Cooper-paired phases. For instance, for QCD with a large number of colors, we know that the so-called chiral density wave Deryagin:1992rw, which is based on fermion-hole pairing (as opposed to fermion-fermion pairing) becomes a strong contender. It has been estimated, however, that at $\mu\sim 1\,{\rm GeV}$ this phase is preferred over color superconductivity only for $N_{c}\gtrsim 1000$ Shuster:1999tn. To find the ground state, we also must compare different color-superconducting phases, which we will be able to do with the result derived in this section.

The free energy density is calculated from Eq. (6). Let us calculate the two terms in this equation separately. As for the gap equation, we shall continue working in the limit of massless quarks and neglect the diagonal components of the quark self-energy for simplicity, $\Sigma^{\pm}\simeq 0$. With ${\rm Tr}\ln=\ln{\rm det}$ and

$${\rm det}=\left(\begin{array}[]{cc}A&B\\ C&D\end{array}\right)={\rm det}(AD-BD^{-1}CD)$$ (59)

for matrices $A,B,C$ and an invertible matrix $D$, we find with Eq. (19)

$$\frac{1}{2}{\rm Tr}\ln\frac{S^{-1}}{T^{2}}=\frac{1}{2}{\rm Tr}\ln\frac{[G_{0}^{+}]^{-1}[G_{0}^{-}]^{-1}-\Phi^{-}G_{0}^{-}\Phi^{+}[G_{0}^{-}]^{-1}}{T^{2}}\,.$$ (60)

We have thus performed the trace over Nambu-Gorkov space, i.e., the trace on the right-hand side is over color, flavor, and Dirac space. We have also included the scale $T^{2}$ in the logarithm, which was not needed in the derivation of the gap equation. Now, with Eqs. (9), (11), and (15) we have

$$[G_{0}^{+}]^{-1}[G_{0}^{-}]^{-1}=\sum_{e}[k_{0}^{2}-(\mu-ek)^{2}]\Lambda_{k}^{-e}\,,$$ (61)

and

$$\Phi^{-}G_{0}^{-}\Phi^{+}[G_{0}^{-}]^{-1}=\Delta^{2}{\cal M}^{2}=\Delta^{2}\sum_{r}\lambda_{r}{\cal P}_{r}\,.$$ (62)

Consequently, we can write

$$\frac{1}{2}{\rm Tr}\ln\frac{S^{-1}}{T^{2}}=\frac{1}{2}{\rm Tr}\ln\sum_{e,r}\Lambda_{k}^{-e}{\cal P}_{r}\frac{k_{0}^{2}-(\epsilon_{k,r}^{e})^{2}}{T^{2}}\,.$$ (63)

The trace of the logarithm of a matrix is the sum over the logarithms of the eigenvalues of the matrix, weighted with their degeneracy. Having written the matrix in terms of projectors, this trace can thus be evaluated easily. With the Dirac trace ${\rm Tr}[\Lambda_{k}^{-e}]=2$ we obtain

$$\frac{1}{2}{\rm Tr}\ln\frac{S^{-1}}{T^{2}} = \frac{T}{V}\sum_{K}\sum_{e,r}{\rm Tr}({\cal P}_{r})\ln\frac{k_{0}^{2}-(\epsilon_{k,r}^{e})^{2}}{T^{2}}$$ (64)

$$= \sum_{e,r}\int\frac{d^{3}\bm{k}}{(2\pi)^{3}}{\rm Tr}({\cal P}_{r})\left[\epsilon_{k,r}^{e}+2T\ln\left(1+e^{-\epsilon_{k,r}^{e}/T}\right)\right]\,,$$

where we have taken the thermodynamic limit and employed the Matsubara sum

$$T\sum_{k_{0}}\ln\frac{k_{0}^{2}-\epsilon_{k}^{2}}{T^{2}}=\epsilon_{k}+2T\ln\left(1+e^{-\epsilon_{k}/T}\right)+{\rm const.}\,,$$ (65)

where the (infinite) terms constant in $k$, $\mu$, and $T$ can be ignored.

For the second term of the free energy density (6) we use Eqs. (8) and (20) to perform the Nambu-Gorkov trace,

$${\rm Tr}(1-S_{0}^{-1}S)={\rm Tr}(2-[G_{0}^{+}]^{-1}G^{+}-[G_{0}^{-}]^{-1}G^{-})\,.$$ (66)

With Eqs. (9) and (26a) we have

$$[G_{0}^{\pm}]^{-1}G^{\pm}=\sum_{e,r}\left(1+\frac{\lambda_{r}\Delta^{2}}{k_{0}^{2}-(\epsilon_{k,r}^{e})^{2}}\right)\Lambda_{k}^{\mp e}{\cal P}_{r}\,.$$ (67)

Therefore,

$$\frac{1}{4}{\rm Tr}(1-S_{0}^{-1}S) = -\frac{T}{V}\sum_{K}\sum_{e,r}{\rm Tr}({\cal P}_{r})\frac{\lambda_{r}\Delta^{2}}{k_{0}^{2}-(\epsilon_{k,r}^{e})^{2}}=\sum_{e,r}\int\frac{d^{3}\bm{k}}{(2\pi)^{3}}{\rm Tr}({\cal P}_{r})\frac{\lambda_{r}\Delta^{2}}{2\epsilon_{k,r}^{e}}\tanh\frac{\epsilon_{k,r}^{e}}{2T}\,,$$ (68)

where we used the Matsubara sum from Eq. (41). Putting both terms (64) and (68) together, we have

$$\Omega = -\sum_{e,r}\int\frac{d^{3}\bm{k}}{(2\pi)^{3}}{\rm Tr}({\cal P}_{r})\left[\epsilon_{k,r}^{e}+2T\ln\left(1+e^{-\epsilon_{k,r}^{e}/T}\right)-\frac{\lambda_{r}\Delta^{2}}{2\epsilon_{k,r}^{e}}\tanh\frac{\epsilon_{k,r}^{e}}{2T}\right]$$ (69)

$$= -\sum_{e,r}\int\frac{d^{3}\bm{k}}{(2\pi)^{3}}{\rm Tr}({\cal P}_{r})\left(\epsilon_{k,r}^{e}-\frac{\lambda_{r}\Delta_{0}^{2}}{2\epsilon_{k,r}^{e}}\right)\,,$$

where, in the second step, we have taken the zero-temperature limit (since $\epsilon_{k,r}^{e}>0$, the logarithm vanishes in this limit). In the absence of a gap, $\Delta_{0}=0$, we find

$$\Omega_{0}=-\frac{N_{c}N_{f}\mu^{4}}{12\pi^{2}}\,,$$ (70)

where we have subtracted the (infinite) vacuum contribution $\Omega(\mu=T=0)$. The result $\Omega_{0}$ is the free energy density of a non-interacting gas of massless quarks because in the present approximation the only effect of interactions is to create a nonzero gap. To compute the contribution of the Cooper pair condensate to the free energy we assume the gap to be nonzero only in the small vicinity of the Fermi surface $k\in[\mu-\delta,\mu+\delta]$ and in this vicinity to assume the constant value $\Delta_{0}$. Then, setting the anti-particle gap to zero, and using

$$\int_{\mu-\delta}^{\mu+\delta}dk\,k^{2}\left(\epsilon_{k,r}^{+}-\frac{\lambda_{r}\Delta_{0}^{2}}{2\epsilon_{k,r}}-|k-\mu|\right)=\frac{\mu^{2}}{2}\left[\lambda_{r}\Delta_{0}^{2}+{\cal O}(\Delta_{0}^{2}/\delta^{2})\right]\,,$$ (71)

we compute

$$\Omega\simeq\Omega_{0}-\frac{\mu^{2}}{4\pi^{2}}\sum_{r}{\rm Tr}({\cal P}_{r})\lambda_{r}\Delta_{0}^{2}\,.$$ (72)

The second term, which is negative, is called condensation energy (density) and shows that the free energy density is indeed reduced by the formation of a Cooper pair condensate. At weak coupling, where $\Delta\ll\mu$, the condensation energy is well approximated by the term of the form $\mu^{2}\Delta_{0}^{2}$. We may use the result to determine the preferred pairing pattern. For CFL, $\sum_{r}{\rm Tr}({\cal P}_{r})\lambda_{r}=8\times 1+1\times 4=12$, while for 2SC $\sum_{r}{\rm Tr}({\cal P}_{r})\lambda_{r}=4\times 1+5\times 0=4$. This number measures how many quasi-particles are gapped and what their gap is in multiples of $\Delta_{0}$. In addition, we need to take into account that $\Delta_{0}$ may vary from phase to phase. We have seen that due to the two-gap structure of CFL, $\Delta_{0}$ is smaller by a factor of $2^{-1/3}\simeq 0.79$ compared to the 2SC phase. This gives a ratio of condensation energies between CFL and 2SC of $3\times 2^{-2/3}\simeq 1.89$. As a consequence, CFL is the favored phase of three-flavor QCD at asymptotically large $\mu$. How far down in density CFL persists and if it (or a variant of it) is still preferred at neutron star densities is unknown. We shall discuss this question on a qualitative level in Sec. III.2.

Without Cooper pairing, the free energy density (or pressure $P=-\Omega$) is of course known beyond the non-interacting limit. More precisely, the pressure is known up to ${\cal O}(\alpha_{s}^{2})$ without Freedman:1976ub and with strange quark mass Kurkela:2009gj, and it is an ongoing effort to compute the complete ${\cal O}(\alpha_{s}^{3})$ contributions Gorda:2021znl; Fernandez:2021jfr; Gorda:2023mkk. It is a challenge to combine these perturbative calculations with the non-perturbative nature of Cooper pairing. Studies in this direction have been made to compute the high-density behavior of the speed of sound Geissel:2024nmx; Geissel:2025vnp.

II.4 Meissner masses and rotated electromagnetism

One of the defining properties of a superconductor is the Meissner effect, i.e., the expulsion of magnetic fields. In field-theoretical terms, the Meissner effect is described by a non-zero magnetic screening mass (the ‘Meissner mass’), whose inverse is the penetration depth, up to which magnetic fields penetrate the superconductor. We can generalize this concept to a color superconductor and compute the magnetic screening masses for the gauge bosons, which, in this case, include the gluons. The screening masses are calculated from the gluon and photon self-energies, the polarization tensors. (For a pedagogical presentation of the calculation in the case of pure electromagnetism, i.e., without gluons, in the same formalism as used here see Ref. Schmitt:2014eka.) At one-loop order the self energies are given by the diagrams in Fig. 5, whose algebraic form is

$$\Pi_{\mu\nu}^{ab}(P)=\frac{1}{2}\frac{T}{V}\sum_{K}{\rm Tr}[\hat{\Gamma}_{a}^{\mu}S(K)\hat{\Gamma}_{b}^{\nu}S(K-P)]\,.$$ (73)

Here, we have introduced the vertices in Nambu-Gorkov space $\hat{\Gamma}_{a}^{\mu}={\rm diag}(\Gamma^{\mu}_{a},\bar{\Gamma}_{a}^{\mu})$, where $\Gamma_{a}^{\mu}=g\gamma^{\mu}T_{a}$ for $a=1,\ldots,8$ and $\Gamma_{a}^{\mu}=e\gamma^{\mu}Q$ for $a=9$, while for the charge-conjugate components $\bar{\Gamma}_{a}^{\mu}=-g\gamma^{\mu}T_{a}^{T}$ for $a=1,\ldots,8$ and $\bar{\Gamma}_{a}^{\mu}=-e\gamma^{\mu}Q$ for $a=9$, where $Q$ is the generator of the electromagnetic gauge group $U(1)_{Q}$ and $e$ is the elementary charge. For the following calculation it is convenient to put the quark flavors in the order $d$, $s$, $u$, such that $Q={\rm diag}(-1/3,-1/3,2/3)=-2/\sqrt{3}\,T_{8}$.

Figure 5: One-loop diagrams for the gluon and photon self-energy in a color superconductor, from which the screening masses are calculated. Solid lines are quark propagators in Nambu-Gorkov space, such that each diagram receives contributions from normal and anomalous propagators, see Fig. 6. Curly external legs are gluons while wavy external legs correspond to photons. In particular, the loop with mixed photon and gluon legs can be nonzero in a color superconductor.

In the combined gluon/photon space, the polarization tensor is a $9\times 9$ matrix and potentially has off-diagonal components that mix gluons with the photon, see Fig. 5. If the full momentum dependence is kept, the polarization tensor can be used to compute the spectral density of the gauge bosons Rischke:2001py; Rischke:2002rz; Malekzadeh:2006ud, as briefly discussed for the HDL approximation in Sec. II.2. While this approximation was sufficient for the calculation of the gap up to a certain order, now we are interested in the effect of Cooper pairing on the gauge bosons themselves. We shall not discuss the full spectral density but are rather interested in the limit $p_{0}=0,p\to 0$. This limit yields the electric (‘Debye’) and magnetic (‘Meissner’) screening masses, which are defined through the longitudinal and transverse components of the self-energy [see Eqs. (29) and (31)],

$$m_{D,ab}^{2} = -\lim_{p\to 0}\Pi^{00}_{ab}(0,\bm{p})\,,$$ (74a)

$$m_{M,ab}^{2} = \frac{1}{2}\lim_{p\to 0}(\delta_{ij}-\hat{p}_{i}\hat{p}_{j})\Pi^{ij}_{ab}(0,\bm{p})\,.$$ (74b)

In the following we will focus on the calculation of the Meissner masses in the CFL phase. We will comment on the Debye masses and other color-superconductors briefly at the end of this section.

Inserting the quark propagator (20) into Eq. (73) and performing the trace over Nambu-Gorkov space yields

$$\Pi^{\mu\nu}_{ab}(P) = \frac{1}{2}\frac{T}{V}\sum_{K}{\rm Tr}[\Gamma_{a}^{\mu}G^{+}(K)\Gamma_{b}^{\nu}G^{+}(Q)+\bar{\Gamma}_{a}^{\mu}G^{-}(K)\bar{\Gamma}_{b}^{\nu}G^{-}(Q)$$ (75)

$$+\Gamma_{a}^{\mu}F^{-}(K)\bar{\Gamma}_{b}^{\nu}F^{+}(Q)+\bar{\Gamma}_{a}^{\mu}F^{+}(K)\Gamma_{b}^{\nu}F^{-}(Q)]\,,$$

where we have abbreviated $Q\equiv K-P$. We see that there are loops formed of normal propagators $G^{\pm}$ and loops formed of anomalous propagators $F^{\pm}$, see Fig. 6.

Figure 6: The gluon (shown here), photon, and mixed polarization tensors receive contributions from normal propagators (left) and anomalous propagators (right), see Eq. (75). The anomalous propagators are represented according to $F^{\pm}=-G_{0}^{\mp}\Phi^{\pm}G^{\pm}$, as in the diagram for the gap equation, Fig. 4. Each diagram also appears with $+$ and $-$ exchanged.

Mixed terms, containing one normal and one anomalous propagator, do not appear. Next, we insert the propagators $G^{\pm}$ and $F^{\pm}$ from Eq. (26). For the spatial components, $\mu=i,\nu=j$, which are needed to compute the Meissner mass, we need the Dirac traces

$${\cal T}^{ij}_{e_{1}e_{2}}(\hat{\bm{k}},\hat{\bm{q}}) \equiv {\rm Tr}[\gamma^{i}\gamma^{0}\Lambda_{k}^{\mp e_{1}}\gamma^{j}\gamma^{0}\Lambda_{q}^{\mp e_{2}}]={\rm Tr}[\gamma^{i}\gamma^{5}\Lambda_{k}^{\pm e_{1}}\gamma^{j}\gamma^{5}\Lambda_{q}^{\mp e_{2}}]$$ (76)

$$= \delta^{ij}(1-e_{1}e_{2}\hat{\bm{k}}\cdot\hat{\bm{q}})+e_{1}e_{2}(\hat{k}^{i}\hat{q}^{j}+\hat{k}^{j}\hat{q}^{i})\,.$$

The color and flavor traces are performed with the help of Eqs. (24) and (25). With ${\rm Tr}[T_{a}T_{b}]=\delta_{ab}/2$ and ${\rm Tr}[T_{a}]=0$ we find for $a,b\leq 8$

$${\cal U}_{ab}^{rs} \equiv g^{2}{\rm Tr}[T_{a}{\cal P}_{r}T_{b}{\cal P}_{s}]=g^{2}{\rm Tr}[T_{a}^{T}{\cal P}_{r}T_{b}^{T}{\cal P}_{s}]=g^{2}\left\{\begin{array}[]{cc}\displaystyle{\frac{7\delta_{ab}}{6}}&\mbox{for}\;\;r=s=1\\[8.61108pt] 0&\mbox{for}\;\;r=s=2\\[8.61108pt] \displaystyle{\frac{\delta_{ab}}{6}}&\mbox{for}\;\;r\neq s\end{array}\right.\,,$$ (80)

and

$${\cal V}_{ab}^{rs} \equiv g^{2}{\rm Tr}[T_{a}{\cal M}{\cal P}_{r}T_{b}^{T}{\cal M}{\cal P}_{s}]=g^{2}{\rm Tr}[T_{a}^{T}{\cal M}{\cal P}_{r}T_{b}{\cal M}{\cal P}_{s}]=g^{2}\left\{\begin{array}[]{cc}0&\mbox{for}\;\;r=s=2\\[8.61108pt] \displaystyle{-\frac{\delta_{ab}}{3}}&\mbox{otherwise}\end{array}\right.\,.$$ (83)

The results for $a=9$ and/or $b=9$ are easily obtained from these expressions: Since $Q=-2/\sqrt{3}\,T_{8}$ we simply multiply these results by $-\frac{2e}{\sqrt{3}g}$ for the traces where one $SU(3)$ generator is replaced by $Q$ and by $\frac{4e^{2}}{9g^{2}}$ if both $SU(3)$ generators are replaced by $Q$. Hence we can write the spatial components of the polarization tensor as

$$\Pi^{ij}_{ab}=\frac{T}{V}\sum_{e_{1},e_{2}}\sum_{r,s}\sum_{K}{\cal T}^{ij}_{e_{1}e_{2}}(\hat{\bm{k}},\hat{\bm{q}})\frac{{\cal U}_{ab}^{rs}(k_{0}q_{0}+e_{1}e_{2}\xi_{k}^{e_{1}}\xi_{q}^{e_{2}})+{\cal V}_{ab}^{rs}\Delta^{2}}{[k_{0}^{2}-(\epsilon_{k,r}^{e_{1}})^{2}][q_{0}^{2}-(\epsilon_{q,s}^{e_{2}})^{2}]}\,,$$ (84)

where we have defined

$$\xi_{k}^{e}\equiv k-e\mu\,.$$ (85)

Next we perform the sum over $k_{0}=-i\omega_{n}$, where $\omega_{n}$ are fermionic Matsubara frequencies. To this end, we make use of the general result

$$T\sum_{k_{0}}\frac{k_{0}q_{0}+\xi_{1}\xi_{2}+\Delta^{2}}{(k_{0}^{2}-\epsilon_{1}^{2})(k_{0}^{2}-\epsilon_{2}^{2})} = \frac{\epsilon_{1}\epsilon_{2}+\xi_{1}\xi_{2}+\Delta^{2}}{2\epsilon_{1}\epsilon_{2}}\left[\frac{f(\epsilon_{1})-f(\epsilon_{2})}{p_{0}-\epsilon_{1}-\epsilon_{2}}-\frac{f(\epsilon_{1})-f(\epsilon_{2})}{p_{0}+\epsilon_{1}+\epsilon_{2}}\right]$$ (86)

$$+\frac{\epsilon_{1}\epsilon_{2}-\xi_{1}\xi_{2}-\Delta^{2}}{2\epsilon_{1}\epsilon_{2}}\left[\frac{1-f(\epsilon_{1})-f(\epsilon_{2})}{p_{0}-\epsilon_{1}-\epsilon_{2}}-\frac{1-f(\epsilon_{1})-f(\epsilon_{2})}{p_{0}+\epsilon_{1}+\epsilon_{2}}\right]$$

$$\to-\frac{\epsilon_{1}\epsilon_{2}-\xi_{1}\xi_{2}-\Delta_{0}^{2}}{2\epsilon_{1}\epsilon_{2}(\epsilon_{1}+\epsilon_{2})}$$

where $q_{0}=k_{0}-p_{0}$ and $p_{0}=-i\omega_{m}$ with the bosonic Matsubara frequencies $\omega_{m}=2\pi mT$, where $f$ is the Fermi-Dirac distribution, and where, in the second step, we have set $p_{0}=0$ and taken the zero-temperature limit, using that $\epsilon_{1},\epsilon_{2}>0$. Moreover, after taking the thermodynamic limit, we need the angular integral

$$\int\frac{d\Omega}{4\pi}{\cal T}^{ij}_{e_{1}e_{2}}(\hat{\bm{k}},\hat{\bm{k}})=\delta^{ij}\left(1-\frac{e_{1}e_{2}}{3}\right)\,.$$ (87)

Putting everything together, we obtain

$$\Pi^{ij}_{ab}(0,\bm{p}\to 0) = -\frac{\delta^{ij}}{3\pi^{2}}\sum_{e}\sum_{r,s}\int_{0}^{\infty}dk\,k^{2}\left\{\frac{{\cal U}_{ab}^{rs}(\epsilon_{k,r}^{e}\epsilon_{k,s}^{e}-\xi_{k}^{e}\xi_{k}^{e})-{\cal V}_{ab}^{rs}\Delta_{0}^{2}}{2\epsilon_{k,r}^{e}\epsilon_{k,s}^{e}(\epsilon_{k,r}^{e}+\epsilon_{k,s}^{e})}+2\frac{{\cal U}_{ab}^{rs}(\epsilon_{k,r}^{e}\epsilon_{k,s}^{-e}-\xi_{k}^{e}\xi^{-e})-{\cal V}_{ab}^{rs}\Delta_{0}^{2}}{2\epsilon_{k,r}^{e}\epsilon_{k,s}^{-e}(\epsilon_{k,r}^{e}+\epsilon_{k,s}^{-e})}\right\}\,.\hskip 14.22636pt$$ (88)

It remains to perform the momentum integral. For now we replace the upper boundary with a finite cutoff $\Lambda$. We first consider the integral over the terms proportional to ${\cal U}_{ab}^{rs}$. We have to distinguish the cases where $r=s$, i.e., where the excitation energies belong to the same quasi-particle branch, and where $r\neq s$, i.e., the terms that mix the two quasi-particle branches with different gaps. For $r=s$ we have

$$\int_{0}^{\Lambda}dk\,k^{2}\left[\frac{(\epsilon_{k,r}^{+})^{2}-(\xi_{k}^{+})^{2}}{2(\epsilon_{k,r}^{+})^{3}}+\frac{(\epsilon_{k,r}^{-})^{2}-(\xi_{k}^{-})^{2}}{2(\epsilon_{k,r}^{-})^{3}}+4\frac{\epsilon_{k,r}^{+}\epsilon_{k,r}^{-}+\xi_{k}^{+}\xi_{k}^{-}}{\epsilon_{k,r}^{+}\epsilon_{k,r}^{-}(\epsilon_{k,r}^{+}+\epsilon_{k,r}^{-})}\right]=2\Lambda^{2}-\mu^{2}+2\Delta_{0}^{2}-3\Delta_{0}^{2}\ln\frac{2\Lambda}{\Delta_{0}}+{\cal O}\left(\frac{1}{\Lambda^{2}}\right).$$ (89)

Although the integrand looks complicated one can perform this integral analytically without approximations. After doing so we have dropped all terms that vanish for $\Lambda\to\infty$. We have produced a divergent term proportional to $\Lambda^{2}$. This is expected and can be removed by subtracting the vacuum contribution. We have also produced a logarithmic divergence that depends on the gap $\Delta_{0}$. If we choose the cutoff to be of the order of $\mu$, say $\Lambda=2\mu$, then we find that this term is exponentially suppressed due to the exponential suppression of the gap. We will thus omit this term in the following. Its appearance is an artifact of the assumption of a momentum-independent gap when we performed the momentum integral. Had we used the full momentum dependence, we would expect a natural cutoff since the gap is peaked at the Fermi surface. As a result, we can approximate the integral (89) by $-\mu^{2}$ since $\mu\gg\Delta_{0}$. With the same arguments we compute the other integrals: for $r\neq s$ we find the same leading order result, such that

$$\sum_{e}\int_{0}^{\infty}dk\,k^{2}\left[\frac{\epsilon_{k,r}^{e}\epsilon_{k,s}^{e}-\xi_{k}^{e}\xi_{k}^{e}}{2\epsilon_{k,r}^{e}\epsilon_{k,s}^{e}(\epsilon_{k,r}^{e}+\epsilon_{k,s}^{e})}+2\frac{\epsilon_{k,r}^{e}\epsilon_{k,s}^{-e}-\xi_{k}^{e}\xi^{-e}}{2\epsilon_{k,r}^{e}\epsilon_{k,s}^{-e}(\epsilon_{k,r}^{e}+\epsilon_{k,s}^{-e})}\right]\simeq-\mu^{2}\,,$$ (90)

independent of $\lambda_{r}$ and $\lambda_{s}$. For the integrals proportional to ${\cal V}_{ab}^{rs}$ (originating from the anomalous propagators) we find

$$\sum_{e}\int_{0}^{\infty}dk\,k^{2}\left[\frac{\Delta_{0}^{2}}{2\epsilon_{k,r}^{e}\epsilon_{k,s}^{e}(\epsilon_{k,r}^{e}+\epsilon_{k,s}^{e})}+2\frac{\Delta_{0}^{2}}{2\epsilon_{k,r}^{e}\epsilon_{k,s}^{-e}(\epsilon_{k,r}^{e}+\epsilon_{k,s}^{-e})}\right]\simeq\left\{\begin{array}[]{cc}\mu^{2}&\mbox{for}\;\;r=s\\[8.61108pt] \displaystyle{\frac{\mu^{2}}{\lambda_{r}-\lambda_{s}}\ln\frac{\lambda_{r}}{\lambda_{s}}}&\mbox{for}\;\;r\neq s\end{array}\right.\,.$$ (93)

Inserting these results into Eq. (88), using the traces (80) and (83) and the definition (74b) we obtain for the Meissner masses of gluons and the photon in the CFL phase

$$m_{M,ab}^{2}=\frac{\mu^{2}}{54\pi^{2}}(21-8\ln 2)\left(\begin{array}[]{cccc}g^{2}&&&\\ &\ddots&&\\ &&g^{2}&\displaystyle{-2/\sqrt{3}\,eg}\\ &&\displaystyle{-2/\sqrt{3}\,eg}&\displaystyle{4/3\,e^{2}}\end{array}\right)\,,$$ (94)

where we have omitted all zeros, and the dots indicate that the diagonal entries for $a=b\leq 8$ are identical. Several observations can be made from this result. We see that Cooper pairing of quarks leads to magnetic screening of gauge bosons, as expected from the analogous effect in ordinary superconductors. In the CFL phase, all 8 gluons acquire a Meissner mass, and the Meissner mass of all gluons is the same. One might wonder why there is no dependence on the gap $\Delta_{0}$. After all, the Meissner mass should vanish in the unpaired phase, where $\Delta_{0}=0$. The reason is that we have taken the limit of vanishing gluon momentum $\bm{p}\to 0$ at a fixed nonzero $\Delta_{0}$. Had we first set $\Delta_{0}=0$ and then computed the screening mass we would have found $m_{M,ab}=0$, as expected. Maybe the most striking feature of the result (94) is the presence of off-diagonal elements. This indicates mixing of the eighth gluon with the photon. By an orthogonal transformation ${\cal O}$ of the gauge fields we can diagonalize the gauge boson propagator, such that $A_{\mu}^{a}\Pi_{ab}^{\mu\nu}A_{\nu}^{b}=\tilde{A}_{\mu}^{a}\tilde{\Pi}_{ab}^{\mu\nu}\tilde{A}_{\nu}^{b}$ with rotated gauge fields $\tilde{A}_{\mu}^{a}={\cal O}_{ab}A_{\mu}^{b}$ and $\tilde{\Pi}_{ab}^{\mu\nu}={\cal O}_{ac}{\Pi}_{cd}^{\mu\nu}{\cal O}_{db}^{T}\propto\delta_{ab}$. In general, this transformation depends on the gauge boson momentum. In the zero-energy, low-momentum case considered here this transformation leaves the gauge fields $A_{\mu}^{1},\ldots,A_{\mu}^{7}$ invariant while the eighth gluon $A_{\mu}^{8}$ and the photon $A_{\mu}$ are rotated,

$$\tilde{A}_{\mu}^{8} = \cos\theta\,A_{\mu}^{8}+\sin\theta\,A_{\mu}\,,$$ (95a)

$$\tilde{A}_{\mu} = -\sin\theta\,A_{\mu}^{8}+\cos\theta\,A_{\mu}\,.$$ (95b)

The mixing angle, which diagonalizes the nontrivial $2\times 2$ sector of the Meissner mass matrix (94) is

$$\cos\theta=\frac{3g^{2}}{3g^{2}+4e^{2}}\,.$$ (96)

After diagonalization, we have $\tilde{m}_{M,88}^{2}=(4e^{2}+3g^{2})/6\,m_{M,88}^{2}$ and $\tilde{m}_{M,99}^{2}=0$, i.e., one combination of the gluon and the photon is massless. This means that there exists a ‘rotated gauge boson’ whose magnetic field can penetrate the CFL phase. In other words, there is a certain charge with respect to which the CFL Cooper pairs are neutral. The generator of the corresponding symmetry is given by a linear combination of the generators of the color and electromagnetic gauge groups, $\tilde{Q}={\bf 1}\otimes Q+\frac{2}{\sqrt{3}}T_{8}\otimes{\bf 1}$, see also the CFL symmetry breaking pattern (1). At strong coupling, $g\gg e$, which is relevant for the physics of neutron stars, the mixing angle (96) becomes very small, $\cos\theta\simeq 1$. This means that the magnetic gauge boson that penetrates the CFL phase is almost identical to the original photon.

The above calculation can also be performed for the electric screening masses and for other color superconductors, see Ref. Schmitt:2003aa for a general discussion and a collection of the results for various phases. In the CFL phase, one finds that the electric screening mass squared is 3 times the magnetic screening mass squared for all $a,b$. In particular, the mixing angle is the same in electric and magnetic sectors. This is different in the 2SC phase, where different mixing angles are found. Moreover, due to the the presence of unpaired quarks, the 2SC screening masses are less symmetric than in the CFL phase: the Debye and Meissner masses of the gluons 1–3 vanish, while the gluons 4–7 have the same nonzero screening, which is different from that of the 8th gluon. As in the CFL phase, there is a residual rotated electromagnetic group, such that a certain rotation of the photon and the 8th gluon has a vanishing Meissner mass.

Rotated electromagnetism plays a role for the calculation of the gap in an external magnetic field. For the CFL phase, the gap equation was solved within an NJL model in Refs. Ferrer:2005vd; Noronha:2007wg, and the QCD gap equation for the 2SC phase with external magnetic field was studied in Ref. Yu:2012jn. It is also interesting to ask if and how the rotated electromagnetism could affect observables of a neutron star with a quark matter core. Within Ginzburg-Landau theory with parameters computed from perturbative QCD, one finds that at strong coupling 2SC and CFL become a type-II superconductor (at $g=3.5$ if at the same time $T_{c}\gtrsim 0.05\mu$ Haber:2017oqb). The corresponding critical magnetic fields for the entrance into a phase of magnetic flux tubes is likely too large to be reached in neutron stars. However, a flux tube lattice is conceivable taking into account the evolution of the star: the magnetic field penetrates the core when the core cools into a color superconductor, and it is unknown on which timescale the magnetic field (more precisely, its rotated version according to the mixing with the gluons) gets expelled or if it stays trapped inside flux tubes. Magnetic flux tubes in CFL can come in different variants, further complicated by the broken $U(1)_{B}$ in the CFL phase, which allows for defects that have magnetic flux and baryon circulation at the same time, i.e., they are superfluid vortices and magnetic flux tubes Iida:2002ev; Iida:2004if; Balachandran:2005ev; Eto:2009kg; Vinci:2012mc; Eto:2013hoa; Alford:2016dco; Haber:2017oqb; Haber:2018tqw. In contrast, 2SC is not a superfluid, its magnetic defects were discussed for instance in Ref. Haber:2017oqb; Alford:2010qf, including magnetic domain walls that exist in the massless limit and turn into multi-winding flux tubes in the presence of a strange quark mass Evans:2020uui. If CFL or 2SC flux tubes are present in a compact star, they may stabilize ‘mountains’ of the star, resulting in detectable continuous gravitational waves Glampedakis:2012qp. Rotated electromagnetism in CFL and in particular the resulting flux tubes also play an important role in the question of a potential quark-hadron continuity. This continuity was originally pointed out based on symmetries and the fermionic spectrum Schafer:1998ef; Alford:1999pa. New aspects arise when flux tubes are taken into account, which is relevant for instance for a potential CFL-hadron interface in a (rotating) neutron star Alford:2018mqj; Chatterjee:2018nxe; Cherman:2018jir; Cherman:2020hbe.

III Color superconductors under neutron star conditions

III.1 Splitting of Fermi momenta

The calculations of the previous subsections were all performed under the premise of weak coupling, treating the coupling constant $g$ as a small parameter. We have briefly discussed extrapolations to large values of $g$ and have mentioned approaches to understand this regime. Besides the problem of strong coupling there is another assumption we have made so far in our discussion of color superconductors. We have assumed all quark masses to vanish. Just like the assumption of the coupling strength to be weak, this assumption is also certain to fail as we move towards lower densities. The reason is that as we go to smaller densities, the relative magnitude of the (strange) quark mass, compared to the chemical potential, increases. This is, firstly, because the chemical potential is smaller at lower densities and, secondly, also because we can expect the (medium-dependent) mass to increase due to strong coupling effects. At moderate densities we expect the strange quark mass to lie somewhere between its current and constituent values, $m_{s}\sim(100-500)\,{\rm MeV}$. This is not negligible compared to the quark chemical potential in the interior of neutron stars $\mu\lesssim 500\,{\rm MeV}$. A difference in quark masses generates a difference in Fermi momenta and thus the usual mechanism of Cooper pairing is challenged. Let us discuss this problem in some more detail, but without going through any full-fledged calculations.

Let us assume that neutron star matter is electrically neutral at every point in space and is in $\beta$-equilibrium, i.e., in equilibrium with respect to the quark analogues of $\beta$-decay and electron capture. (Ignoring mixed phases, where local neutrality can be violated Glendenning:1992vb; Heiselberg:1992dx; Schmitt:2020tac, and ignoring fast dynamics as in neutron star mergers, where matter may be out of chemical equilibrium Hammond:2021vtv; Chabanov:2023abq.) This imposes a constraint on quark and electron densities,

$$\frac{2}{3}n_{u}-\frac{1}{3}n_{d}-\frac{1}{3}n_{s}-n_{e}=0\,,$$ (97)

and on the chemical potentials

$$\mu_{u}+\mu_{e}=\mu_{d}\,,\qquad\mu_{u}+\mu_{e}=\mu_{s}\,,$$ (98)

where we have neglected the neutrino chemical potential, assuming that the neutrino mean free path is of the order of or larger than the size of the star. We also have to impose color neutrality, although we will largely ignore this point in the following, mostly qualitative, discussion. Massless three-flavor quark matter is particularly symmetric: equal quark densities and chemical potentials fulfill these conditions, and no electrons are needed. The reason is that the electric charges of the three lightest quarks happen to add up to zero. Therefore, in the previous subsections it was justified to consider a single, common quark chemical potential $\mu$ for all quark species. If we now ‘switch on’ a mass for the strange quark, its Fermi momentum is reduced to $k_{F,s}=\sqrt{\mu_{s}^{2}-m_{s}^{2}}$; it has become energetically more costly to populate the system with strange quarks. Therefore, since $n_{s}\propto k_{F,s}^{3}$ (at $T=0$), the number density of strange quarks is reduced as well. As a consequence, and because the up and down quarks remain approximately massless, we have disrupted the symmetric situation and the quark and electron densities – under the conditions (97) and (98) – will adjust themselves in a less symmetric way. If all interactions are neglected, one computes that, to lowest order in the strange quark mass,

$$\mu_{e}\simeq\frac{m_{s}^{2}}{4\mu}\,,$$ (99)

and

$$k_{F,u}\simeq\mu-\frac{m_{s}^{2}}{6\mu}\,,\qquad k_{F,d}\simeq\mu+\frac{m_{s}^{2}}{12\mu}\,,\qquad k_{F,s}\simeq\mu-\frac{5m_{s}^{2}}{12\mu}\,.$$ (100)

We see that the Fermi momenta of the quarks are split by an equal amount of

$$\delta k_{F}\simeq\frac{m_{s}^{2}}{4\mu}\,,$$ (101)

such that $k_{F,s}<k_{F,u}<k_{F,d}$. The lack of negative charge due to the smaller number of strange quarks is compensated by an increased Fermi momentum of the down quark and by the presence of electrons. (Note that ‘the down goes up and the up goes down’.) The split $\delta k_{F}$ becomes larger as we go down in density since $m_{s}$ becomes larger and $\mu$ becomes smaller.

III.2 Cooper pairing beyond CFL

How does a split of Fermi momenta affect Cooper pairing? This is a very general question, which, in its simplest form, is: consider a (zero-temperature) system of two fermion species with attractive interaction and chemical potentials $\mu_{1}$ and $\mu_{2}$, such that the mismatch in chemical potentials is $\delta\mu=(\mu_{1}-\mu_{2})/2$ and the average chemical potential is $\bar{\mu}=(\mu_{1}+\mu_{2})/2$. Is ‘usual’ Cooper pairing possible at a nonzero $\delta\mu$ and if yes at what value of $\delta\mu$ does it break down? The answer is that it is still energetically favorable to form Cooper pairs if the mismatch $\delta\mu$ is sufficiently small compared to the gap $\Delta$. In this case, the pairing process can be viewed as follows. Starting from the two different Fermi spheres, first create a fictitious state in which both fermion species are filled up to a common Fermi momentum. This costs free energy. Then pair. If the energy gain from pairing is larger than the cost from creating a common Fermi momentum, pairing will happen. In this case, the quasi-fermion dispersion relations are modified by the mismatch to $\epsilon_{k}^{e}\pm\delta\mu$ with $\epsilon_{k}^{e}$ from Eq. (27). In particular, if $\delta\mu>\Delta$, one excitation branch becomes gapless despite the presence of a nonzero $\Delta$. However, one finds that Cooper pairing breaks down before gapless modes appear. The critical value for the mismatch beyond which the usual pairing is no longer favored over the unpaired phase turns out to be $\delta\mu=\Delta/\sqrt{2}$. This is called the Chandrasekhar-Clogston limit Chandrasekhar:1962; Clogston:1962. For a pedagogical discussion of the calculation leading to this result see for instance Ref. Schmitt:2014eka.

In quark matter, the situation is more complicated because we have 9 fermion species, not just two, and we have the conditions of neutrality and $\beta$-equilibrium. In QCD, the split between the Fermi surfaces as well as the gap $\Delta$ are all functions of the single thermodynamic parameter $\mu$ (at $T=0$). Nevertheless, the general principle is the same: a mismatch in Fermi momenta tends to break Cooper pairs. And, it is not possible to avoid the mismatch and maintain ‘conventional pairing’ by choosing a different pairing pattern Rajagopal:2005dg. It turns out that CFL becomes energetically disfavored compared to the unpaired phase at $m_{s}^{2}/\mu=4\Delta$ Alford:2002kj. However, one finds that gapless excitations already set in at the smaller mismatch $m_{s}^{2}/\mu=2\Delta$. (Gaplessness first occurs in the sector where $bd$ quarks pair with $gs$ quarks. Due to the presence of color chemical potentials the mismatch for this sector is $\delta\mu=m_{s}^{2}/(2\mu)$, not $m_{s}^{2}/(4\mu)$, as the result of unpaired quark matter (101) would suggest.) Although, by comparing the free energy with the unpaired phase, ‘gapless CFL’ appears to be stable, it suffers an instability with respect to the formation of anisotropic or even inhomogeneous phases. This instability manifests itself in a negative Meissner mass squared, i.e., if the calculation of Sec. II.4 is repeated in the presence of a strange quark mass, one finds unphysical Meissner masses as soon as at least one excitation branch becomes gapless Casalbuoni:2004tb; Fukushima:2005cm. This was first pointed out in the context of the 2SC phase Huang:2004bg.

These results give us some idea about whether we can expect CFL in the interior of neutron stars, but due to the lack of knowledge of the medium-dependent $m_{s}$ and the color-superconducting gap $\Delta$ at moderate densities we do not have a conclusive answer to this question. What if CFL breaks down before the transition to nuclear matter is reached? In this case, rather than being completely unpaired, quark matter is expected to be in a more exotic paired phase. Almost all of these phases are ‘partially paired’ phases, where some, but not all, fermion species remain unpaired and/or where certain fermion species only form pairs in certain directions in momentum space. In either case, these phases are interesting for the phenomenology of quark matter in neutron stars because they differ in their transport properties from the completely gapped CFL phase. Let us discuss these phases a bit more systematically.

If the parameter $m_{s}^{2}/(\mu\Delta)$ is small, one can discuss the fate of CFL in a controlled manner within weak coupling. It turns out that the next phase down in density then is the so-called CFL-$K^{0}$ phase, where a condensate of neutral kaons forms on top of the fermionic Cooper pair condensate Bedaque:2001je; Buballa:2004sx; Forbes:2004ww; Alford:2007qa. To understand this phase we recall that CFL breaks the approximate chiral symmetry of QCD, such that there is an octet of pseudo-Goldstone modes. The quantum numbers of these modes are the same as the ones of the usual meson octet in low-density QCD. This is due to the identical symmetry breaking pattern regarding the global chiral symmetry. Therefore, there are CFL versions of the pions, kaons, and the eta, which can be described by the chiral field

$$\Sigma=\phi_{L}^{\dagger}\phi_{R}\in SU(3)\,,$$ (102)

where $\phi_{L}$ and $\phi_{R}$ are $3\times 3$ matrices for the left-handed and right-handed sectors. This is to be compared to Eq. (17), where we assumed $\phi_{L}=\phi_{R}\equiv\phi$. In the ‘pure’ CFL phase, therefore, $\Sigma=1$. Each $\phi$ in Eq. (102) represents diquarks, and thus the CFL mesons are four-fermion objects, made of two fermions and two holes. Therefore, even though their quantum numbers are the same as for the usual mesons, their properties, in particular their masses, are different. For instance, the CFL neutral kaon can be viewed as an object $K^{0}\sim\bar{u}\bar{s}du$. This can be obtained from the ‘usual’ kaon $K^{0}\sim\bar{s}d$ by the replacements $s\to\bar{u}\bar{d}$ and $d\to\bar{u}\bar{s}$. These replacements reflect the anti-triplet, which the flavor index of $\phi_{R}$ and $\phi_{L}$ represents (the second index is a color index). As a consequence, the meson masses are ordered differently: while the quark flavors $(u,d,s)$ have mass ordering $m_{u}<m_{d}<m_{s}$, the corresponding anti-triplet $(\bar{d}\bar{s},\bar{u}\bar{s},\bar{u}\bar{d})$ has the opposite ordering $m_{d}m_{s}>m_{u}m_{s}>m_{u}m_{d}$. Therefore, in the QCD vacuum, $m_{\pi^{0}}<m_{K^{0}}$ since $m_{\pi^{0}}\propto m_{u}+m_{d}$ and $m_{K^{0}}\propto m_{s}+m_{d}$, while in CFL $m_{K^{0}}<m_{\pi^{0}}$ since $m_{K^{0}}\propto m_{u}m_{d}+m_{u}m_{s}$ while $m_{\pi^{0}}\propto m_{d}m_{s}+m_{u}m_{s}$. Indeed, the neutral kaon is the lightest CFL meson. Masses and other properties of CFL mesons can be studied within a chiral effective theory with Lagrangian Son:1999cm; Bedaque:2001je

$${\cal L}=\frac{f_{\pi}^{2}}{4}{\rm Tr}[\nabla_{0}\Sigma\nabla_{0}\Sigma^{\dagger}-v_{\pi}^{2}\partial_{i}\Sigma\partial_{i}\Sigma^{\dagger}]+\frac{af_{\pi}^{2}}{2}{\rm det}\,M\,{\rm Tr}[M^{-1}(\Sigma+\Sigma^{\dagger})]\,,$$ (103)

with $v_{\pi}=1/\sqrt{3}$, the covariant derivative $\nabla_{0}\Sigma\equiv\partial_{0}\Sigma+i[A,\Sigma]$, which includes effective chemical potentials of the mesons through $A\equiv-\frac{M^{2}}{2\mu}$, and the quark mass matrix $M={\rm diag}(m_{u},m_{d},m_{s})$. The constants $f_{\pi}$ and $a$ can be computed at weak coupling,

$$f_{\pi}^{2}=\frac{21-8\ln 2}{18}\frac{\mu^{2}}{2\pi^{2}}\,,\qquad a=\frac{3\Delta^{2}}{\pi^{2}f_{\pi}^{2}}\,.$$ (104)

Since the structure of this effective theory is determined entirely by symmetries, it gives us a rigorous low-energy description of the CFL phase. Notice that all quarks are gapped in ‘pure’ CFL and thus its transport properties are entirely determined by the (pseudo-)Goldstone modes. Thus, if fully gapped CFL survives in neutron star conditions, the chiral Lagrangian (103) provides us with a rigorous tool to study quark matter at moderate densities, up to the uncertainties in the parameters $f_{\pi}$ and $a$. The Lagrangian (103) can be used to demonstrate that a condensate of neutral kaons becomes preferred if the effective kaon chemical potential $\mu_{K^{0}}$ becomes larger than the kaon mass $m_{K^{0}}$. With

$$\mu_{K^{0}}=\frac{m_{s}^{2}-m_{d}^{2}}{2\mu}\,,\qquad m_{K^{0}}^{2}=am_{u}(m_{s}+m_{d})\,,$$ (105)

and the weak-coupling results (104) for $f_{\pi}$ and $a$ we find that the critical strange mass for kaon condensation scales as $m_{u}^{1/3}\Delta^{2/3}$. In this kaon-condensed CFL phase, all fermion excitations remain fully gapped. As we go further down in density, however, ungapped fermions appear. They do so first in certain regions in momentum space in the so-called curCFL-$K^{0}$ phase Schafer:2005ym; Kryjevski:2008zz. In this phase, a meson supercurrent appears, which is balanced by a counter-propagating current of ungapped fermions. This phase thus spontaneously breaks rotational invariance. More counter-propagating currents in more directions may appear as we go further down in density. The resulting phases break translational invariance and exhibit periodically varying gaps Alford:2000ze; Casalbuoni:2001gt; Giannakis:2002jh; Casalbuoni:2005zp; Carignano:2017meb. These crystalline phases are also called Larkin-Ovchinnikov-Fulde-Ferrell (LOFF) phases FuldeFerrell; LarkinOvchinnikov, and different lattice structures have been studied to find the preferred phase Rajagopal:2006ig. These calculations are based on Ginzburg-Landau studies or phenomenological models, and thus we are starting to enter theoretically uncontrolled territory in our journey down in density. Nevertheless, the possible existence of crystalline quark phases is of interest for neutron star physics, for example in the context of pulsar glitches Mannarelli:2007bs, continuous gravitational wave emission Knippel:2009st, or the tidal deformability Lau:2017qtz; Lau:2018mae.

Besides breaking CFL Cooper pairs in some directions in momentum space, the system might also react by leaving strange quarks completely unpaired, which is the 2SC phase. We have seen that at asymptotically large densities the 2SC phase has larger free energy than CFL and thus is not preferred. This remains true if a small strange quark mass is switched on. However, as the coupling and/or the mismatch in chemical potentials is increased, a transition from CFL to 2SC becomes possible. For instance in NJL calculations with strong diquark coupling the 2SC phase tends to occupy a large region of the phase diagram between CFL and nuclear matter Ruester:2005jc; Gholami:2024diy. Since 2SC is also based on cross-flavor pairing, a split between $u$ and $d$ Fermi surfaces tends to break 2SC for the same reasons as discussed for CFL, and a 2SC version of the LOFF phase becomes a candidate phase.

If the mismatch between Fermi surfaces is too large, be it in CFL or 2SC, the system may resort to single-flavor pairing Schafer:2000tw; Schmitt:2004hg; Schmitt:2004et; Aguilera:2005tg; Schmitt:2005wg; Alford:2005yy; Feng:2007bg; Brauner:2008ma; Pang:2010wk; Fujimoto:2019sxg; Fujimoto:2021bes; Sogabe:2024yfl. Due to the anti-symmetry of the Cooper pair wave function [see discussion below Eq. (16)], this requires pairing in the spin-1 channel. As we have seen above, the weak-coupling gap in the spin-1 channel is about 2-3 orders of magnitude smaller than the spin-0 gap. Therefore, the gain in condensation energy is much smaller as well, and spin-1 pairing is only expected to occur if cross-flavor pairing becomes impossible due to a large mismatch in Fermi momenta or as a “residual” pairing of fermionic modes that are otherwise left unpaired, for instance strange quarks Alford:2005yy or blue down quarks Fujimoto:2019sxg in the 2SC phase. Spin-one color superconductivity allows for a multitude of possible pairing patterns in itself. The symmetries involved are not unlike the ones in superfluid ³He, and thus many candidate phases have their analogues in ³He. The energetically preferred phase at weak coupling is the color-spin locked phase Schafer:2000tw; Schmitt:2004et, analogous to the B-phase in ³He Vollhardt1990. This is the only spin-1 color superconductor in which the gap function is isotropic. Moreover, all Meissner masses, including that of the photon are nonzero Schmitt:2003xq. Other spin-1 phases have anisotropic gaps and may have nodal points or lines at the Fermi surface, such as the polar and A phases Schmitt:2004et.

IV Bulk viscosity of color superconductors

IV.1 Transport in neutron stars

Transport plays a crucial role for the understanding of color-superconducting phases. While the equation of state is only mildly sensitive to Cooper pairing, transport properties can change dramatically if matter becomes superconducting or superfluid. All particles in the system contribute to thermodynamic properties such as the pressure, even if the particles are excluded from any scattering processes due to Pauli blocking. Therefore, if a small fraction of particles in a small vicinity of the Fermi surface undergoes Cooper pairing (having in mind a weakly coupled and sufficiently cold system) this does not significantly change the bulk properties such as the pressure. Transport, however, is sensitive to a change in properties of the particles at the Fermi surface, because these are the particles that dominate transport to begin with, through scattering processes or other reactions such as flavor-changing processes. (The stronger the coupling, the less accurate this distinction between bulk physics and the physics at the Fermi surface becomes, and at very strong coupling the whole picture of transport through quasi-particles becomes questionable.)

Here, ‘transport’ describes the transfer of conserved quantities from one spatial region to another. This occurs if the system is brought out of equilibrium and attempts to re-equilibrate. For instance, a temperature gradient will induce a transfer of energy from the hotter to the colder region, or a gradient in chemical potentials will induce processes that change the chemical composition of matter. Calculating ‘transport properties’ therefore amounts to determine how efficiently, i.e., on which time scale, the system is able to re-equilibrate.

There are many observables of neutron stars that are sensitive to transport properties and thus are able to discriminate between paired and unpaired phases or between different (color-)superconducting phases. Examples are neutrino emissivity and the resulting cooling behavior Potekhin:2015qsa; Marino:2024gpm, damping of oscillations and consequences for the rotational frequency of the star Andersson:1997xt; Kraav:2024cus, dissipation in neutron star mergers Alford:2017rxf, pulsar glitches Antonopoulou:2022rpq, and the evolution of the magnetic field Pons:2019zyc.

In a neutron star, there are many particle species and many possible scattering and reaction processes that can potentially contribute to transport, see Ref. Schmitt:2017efp for a review that takes into account all layers of the star, from the crust to a potential quark matter core. Typically, if present, fermionic contributions (electrons, nucleons, quarks, $\ldots$) dominate over bosonic ones (kaons, pions, superfluid phonons, $\ldots$), at least at sufficiently small temperatures, where, due to the presence of the Fermi surface, fermionic contributions are usually enhanced by powers of $\mu/T$ compared to bosonic ones. Fermions, however, may become ‘unavailable’ for low-energy transport through Cooper pairing if the temperature is much smaller than the energy gap induced in the fermionic spectrum due to pairing. In most forms of neutron star matter, not all fermionic degrees of freedom are gapped. For instance in nuclear matter, even if protons and neutrons are paired, unpaired electrons are present and will dominate most transport properties. In quark matter, CFL is special because all quarks are gapped and no electrons are expected at zero temperature. In this case, transport is dominated by mesons and the superfluid phonon Jaikumar:2002vg; Shovkovy:2002sg; Alford:2007rw; Manuel:2007pz; Alford:2008pb; Mannarelli:2009ia; Alford:2009jm; Anglani:2011cw; Bierkandt:2011zp. In contrast, in the 2SC phase, there are ungapped quarks and a sizable fraction of electrons, both of which are relevant for the calculation of transport properties Jaikumar:2005hy; Alford:2006gy; Schmitt:2010pn; Alford:2014doa; Alford:2025tbp; Alford:2025jtm. This is similar in all color-superconducting phases other than CFL, for instance in spin-one color superconductors Schmitt:2005wg; Sad:2006egl; Wang:2009if; Wang:2010ydb; Berdermann:2016mwt or in the LOFF phase Anglani:2006br; Sedrakian:2013pva; Sarkar:2016gib.

We shall use bulk viscosity of the 2SC phase from the processes $u+d\leftrightarrow u+s$ as an example for a transport coefficient of a color-superconducting phase. Going through the details of the calculation will show the complications of calculating transport coefficients in a superconducting phase, and also will allow us to derive a simple analytical result, at least in certain limits, for unpaired quark matter. Bulk viscosity is different from shear viscosity, thermal conductivity, and electric conductivity in the sense that its dominant contribution originates from flavor-changing electroweak processes. These processes are most efficient in creating dissipation because they can occur on the same time scales as given by the oscillations of the star, or by the merger of two neutron stars. These time scales are of the order of milliseconds or larger and thus always much larger than those of scattering processes of the strong interaction. Therefore (and in contrast to heavy-ion collisions) bulk viscosity of neutron star matter is dominated by processes of the electroweak interaction, not by QCD scattering. Strong interaction effects may nevertheless provide corrections to the electroweak rate, for instance through so-called non-Fermi liquid effects, which alter the dispersion relations of (unpaired) quarks Schafer:2004zf; Schwenzer:2012ga. In the following calculation we shall ignore these corrections for simplicity. To make our calculation self-contained, we start with the derivation of the bulk viscosity.

IV.2 Derivation of bulk viscosity

Bulk viscosity is responsible for dissipation due to compression and expansion of a fluid. We denote the angular frequency of the external volume oscillation by $\omega$, such that we can write the volume as $V(t)=V_{0}[1+\delta v(t)]$, with a dimensionless volume perturbation

$$\delta v(t)=\delta v_{0}\cos\omega t\,,$$ (106)

which we assume to be small, $\delta v_{0}\ll 1$. Denoting averaging over one oscillation period $\tau=2\pi/\omega$ by angular brackets, we can write the dissipated energy density, on the one hand, as

$$\langle\dot{{\cal E}}\rangle_{\tau}=-\frac{\zeta}{\tau}\int_{0}^{\tau}dt\,(\nabla\cdot{\bm{v}})^{2}\simeq-\frac{\zeta\omega^{2}\delta v_{0}^{2}}{2}\,,$$ (107)

where $\zeta$ is the bulk viscosity. Here we have used the continuity equation at zero three-velocity ${\bm{v}}=0$ (but nonzero derivative of $\bm{v}$) to relate the divergence of the velocity field to the change in the total particle number density, $n\nabla\cdot\bm{v}=-\partial_{t}n$. With $n=N/V$ and holding the particle number $N$ fixed, this is then rewritten in terms of the change in volume. On the other hand, the dissipated energy density can be expressed in terms of the mechanical work done by the induced pressure oscillations,

$$\langle\dot{{\cal E}}\rangle_{\tau}=\frac{1}{\tau}\int_{0}^{\tau}dt\,P(t)\frac{d\delta v}{dt}\,.$$ (108)

We write the pressure of three-flavor quark matter as

$$P(t)=P_{0}+\frac{\partial P}{\partial V}V_{0}\delta v(t)+\sum_{x=u,d,s}\frac{\partial P}{\partial n_{x}}\delta n_{x}(t)\,,$$ (109)

where the derivatives on the right-hand side are evaluated in equilibrium. The oscillation in the pressure is in general out of phase compared to the volume oscillation because of the microscopic re-equilibration processes which induce changes in the number densities of the particle species $\delta n_{x}$. From Eqs. (107) and (108) we obtain the bulk viscosity in the form

$$\zeta=-\frac{2}{\omega^{2}\delta v_{0}^{2}}\frac{1}{\tau}\int_{0}^{\tau}dt\,P(t)\frac{d\delta v}{dt}\,.$$ (110)

To calculate this expression, we assume that the dominant electroweak re-equilibration processes are the non-leptonic reactions

$$u+d \leftrightarrow u+s\,.$$ (111)

There are other processes that potentially contribute to the bulk viscosity, such as the leptonic processes (‘direct Urca’ processes)

$$u+e \to d+\nu_{e}\,,\qquad d\to u+e+\bar{\nu}_{e}$$ (112a)

$$u+e \to s+\nu_{e}\,,\qquad s\to u+e+\bar{\nu}_{e}\,.$$ (112b)

Their contribution becomes relevant for temperatures larger than typically found in isolated neutron stars, see Ref. Alford:2025tbp and below discussion of Fig. 9. (See Refs. Iwamoto:1980eb; Iwamoto:1982zz for the calculation of the rates of these processes in unpaired quark matter.)

In chemical equilibrium, the two reactions (111) do not change the various densities because they occur with the same rate, and according to the principle of detailed balance in this situation the sum of the chemical potentials of the ingoing particles equals the sum of the chemical potentials of the outgoing particles. Therefore, defining

$$\delta\mu\equiv\mu_{s}-\mu_{d}\,,$$ (113)

chemical equilibrium implies $\delta\mu=0$. A deviation from chemical equilibrium occurs if the equality of chemical potentials is disrupted, $\delta\mu\neq 0$. The situation considered here is particularly simple because there is a single process (111) and a single $\delta\mu$. In general, there can be multiple processes related to the same $\delta\mu$, or there can be multiple processes related to multiple $\delta\mu$’s. For instance, if the leptonic processes (112) were taken into account there would be two relevant differences in chemical potentials, $\delta\mu_{1}=\mu_{s}-\mu_{d}$ and $\delta\mu_{2}=\mu_{d}-\mu_{u}-\mu_{e}$ (neglecting the neutrino chemical potential since neutrinos and anti-neutrinos leave the system once they are created, at least if the temperature is sufficiently low). The resulting expression for the bulk viscosity then becomes significantly more cumbersome, see for instance Ref. Sad:2007afd or appendix A of Ref. Alford:2006gy. Here we proceed with a single process and a single $\delta\mu$, in which case we can write the change in quark number densities due to the electroweak processes (111) as

$$\delta n_{d}(t)=-\delta n_{s}(t)=\int_{0}^{t}dt^{\prime}\,\Gamma[\delta\mu(t^{\prime})]\simeq\lambda\int_{0}^{t}dt^{\prime}\delta\mu(t^{\prime})\,,\qquad\delta n_{u}(t)=0\,.$$ (114)

Here, $\Gamma[\delta\mu(t)]$ is the net number of $d$ quarks produced per unit time and volume in the processes $u+s\leftrightarrow u+d$. Calculating $\Gamma$ is the main task in calculating the bulk viscosity and will be explained in the next subsection. We have linearized the result for small $\delta\mu(t)$, such that now all microscopic input is encoded in $\lambda$. According to our definition of $\delta\mu$, a net production of $d$ quarks sets in for $\delta\mu>0$, from which we conclude $\lambda>0$.

The difference in chemical potentials $\delta\mu$ oscillates due to the volume oscillation and due to the electroweak reactions,

$$\frac{d\delta\mu}{dt} = \frac{\partial\delta\mu}{\partial V}\frac{dV}{dt}+\sum_{x=d,s}\frac{\partial\delta\mu}{\partial n_{x}}\frac{dn_{x}}{dt}\,.$$ (115)

Using Eqs. (113), (114), the thermodynamic equilibrium relation

$$\frac{\partial P}{\partial n_{x}}=-V_{0}\frac{\partial\mu_{x}}{\partial V}\,,$$ (116)

and abbreviating

$$B\equiv\frac{\partial P}{\partial n_{s}}-\frac{\partial P}{\partial n_{d}}\,,\qquad C\equiv\frac{\partial\mu_{s}}{\partial n_{s}}+\frac{\partial\mu_{d}}{\partial n_{d}}-\frac{\partial\mu_{d}}{\partial n_{s}}-\frac{\partial\mu_{s}}{\partial n_{d}}\,,$$ (117)

we obtain

$$\frac{d\delta\mu}{dt} = -B\frac{d\delta v}{dt}-\lambda C\delta\mu(t)\,.$$ (118)

Note that $B$ and $C$ are thermodynamic functions evaluated in equilibrium, i.e., they do not depend on the electroweak reaction rate. The pressure (109) can now be expressed with the help of $B$,

$$P(t)=P_{0}+\frac{\partial P}{\partial V}V_{0}\delta v(t)-B\delta n_{d}(t)\,,$$ (119)

which, with Eq. (108), yields the dissipated energy density in the form

$$\langle\dot{{\cal E}}\rangle_{\tau}=\left\langle P(t)\frac{d\delta v}{dt}\right\rangle=\lambda B\langle\delta\mu(t)\delta v(t)\rangle\,,$$ (120)

where partial integration and Eq. (114) have been employed. Notice that there is only dissipation if $B$ is nonzero. The reason is that $B$ tells us to what extent the system is brought out of equilibrium by a volume change. For instance, for non-interacting quark matter at zero temperature and neglecting the light quark masses,

$$B=\left(n_{s}\frac{\partial\mu_{s}}{\partial n_{s}}-n_{d}\frac{\partial\mu_{d}}{\partial n_{d}}\right)_{\delta\mu=0}=-\frac{m_{s}^{2}}{3\mu_{d}}\,.$$ (121)

In other words, if strangeness does not introduce an asymmetry in the thermodynamic relations in equilibrium, there will be no dissipation and thus no bulk viscosity due to a strangeness changing process.

In general, $\delta\mu(t)$ oscillates out of phase with the volume $\delta v(t)$, and thus we make the ansatz $\delta\mu(t)={\rm Re}[\delta\mu_{0}e^{i\omega t}]$, with a complex amplitude $\delta\mu_{0}$. The differential equation (115) then yields algebraic equations for real and imaginary parts of $\delta\mu_{0}$, which can easily be solved,

$${\rm Re}\,\delta\mu_{0}=-\frac{B\omega^{2}\delta v_{0}}{\omega^{2}+(C\lambda)^{2}}\,,\qquad{\rm Im}\,\delta\mu_{0}=-\frac{BC\lambda\omega\delta v_{0}}{\omega^{2}+(C\lambda)^{2}}\,.$$ (122)

In order to obtain this simple analytical result it was crucial to use the linear approximation of the rate (114). As we shall see, higher powers in $\delta\mu$ do appear in the rate in general. The linear approximation is valid in the so-called ‘sub-thermal’ regime $\delta\mu\ll T$. In the ‘supra-thermal’ regime, however, the deviation from chemical equilibrium becomes of the order of the temperature and larger PhysRevD.46.3290; Alford:2010gw.

With the help of Eq. (122) we can compute $\delta n_{d}$ from the integral in Eq. (114), insert the result into the pressure (119), and the result into the expression for the bulk viscosity (110). This yields

$$\zeta(\omega)=\frac{\lambda B^{2}}{(\lambda C)^{2}+\omega^{2}}\,.$$ (123)

For given values of the external frequency $\omega$ and the equilibrium functions $B$ and $C$, the bulk viscosity has a maximum at $\lambda=\omega/C$. In this sense, bulk viscosity is a resonance phenomenon. In particular, the bulk viscosity relevant for the physics of neutron stars is frequency dependent. The zero-frequency limit is not sufficient because there are re-equilibration processes that occur on the same time scale as the external oscillation. As we shall see, the linearized rate $\lambda$ increases monotonically with temperature (while for $B$ and $C$ we can, to a good approximation, use the zero-temperature limit). This implies that the viscosity for a given external frequency $\omega$ has a maximum at a certain temperature. It turns out that this temperature is of the order of 1 MeV or even larger. Thus, bulk viscosity, in contrast to shear viscosity, is particularly important for young neutron stars and potentially for neutron star mergers. The purpose of the following subsection is to calculate $\lambda$.

IV.3 Calculating the rate of $u+d\leftrightarrow u+s$ in 2SC quark matter

IV.3.1 Setup

Our starting point is the kinetic equation, which arises from the more general Kadanoff-Baym equation Kadanoff via a gradient expansion,

$$i\frac{\partial}{\partial t}{\rm Tr}[\gamma_{0}S^{<}(P_{1})]=-{\rm Tr}[S^{>}(P_{1})\Sigma^{<}(P_{1})-\Sigma^{>}(P_{1})S^{<}(P_{1})]\,,$$ (124)

where $S^{<}$ and $S^{>}$ are the lesser and greater quark propagators in momentum space, $\Sigma^{<}$ and $\Sigma^{>}$ are the lesser and greater quark self-energies, $P_{1}$ is the $d$ quark four-momentum, and the trace is taken over color, flavor, Dirac, and Nambu-Gorkov space. The greater and lesser propagators are an important ingredient of the real-time formalism, see for instance Ref. lebellac. By inserting the equilibrium propagators, the kinetic equation will give us the net number of $d$ quarks produced per unit time and volume $\Gamma$ through the process $u+d\leftrightarrow u+s$. This process enters the equation through the appropriate choice of the self energy. The left-hand side of the kinetic equation essentially is $\Gamma$ (after integrating over momentum and up to a numerical factor), and our task is to compute the right-hand side, which gives the gain and loss terms due to the processes $u+s\to u+d$ and $u+d\to u+s$. Starting from the kinetic equation rather than writing down the rate more directly from Fermi’s Golden Rule is particularly useful for superconducting phases since it naturally yields all sub-processes that are specific to Cooper-paired systems and whose form and even existence are not obvious.

The self-energies $\Sigma^{<,>}$ are given by the diagram shown in Fig. 7.

Figure 7: Left: Quark self-energy needed for the calculation of the rate $u+d\leftrightarrow u+s$. Right: The diagram is duplicated to show our notation for the four-momenta of the various propagators. Due to energy-momentum conservation, $Q=P_{1}-P_{4}=P_{3}-P_{2}$.

Cutting this diagram in half yields the processes we are interested in. Since the $W$-boson mass $M_{W}$ is much larger than any other energy involved, the $W$-boson lines can be reduced to pointlike interactions, and the self-energies are

$$\Sigma^{<,>}(P_{1})=\frac{i}{M_{W}^{4}}\int_{P_{4}}\,\Gamma_{ud,-}^{\mu}\,S^{<,>}(P_{4})\,\Gamma_{ud,+}^{\nu}\,\Pi_{\mu\nu}^{>,<}(Q)\,,$$ (125)

where we have abbreviated the four-momentum integration by

$$\int_{P}\equiv\int\frac{d^{4}P}{(2\pi)^{4}}\,,$$ (126)

and where $\Gamma_{ud,\pm}^{\mu}$ are the vertices of the subprocesses $d\leftrightarrow u+W^{-}$. The $W$-boson polarization tensors are

$$\Pi_{\mu\nu}^{<,>}(Q)=-i\int_{P_{2}}{\rm Tr}[\Gamma_{us,+}^{\mu}\,S^{>,<}(P_{3})\,\Gamma_{us,-}^{\nu}\,S^{<,>}(P_{2})]\,,$$ (127)

where $\Gamma_{us,\pm}^{\mu}$ are the vertices of the subprocesses $u+W^{-}\leftrightarrow s$. In the diagram, the $W$-boson polarization tensor corresponds to the quark loop of $u$ and $s$ quarks. This is analogous to the photon and gluon polarization tensors used in the calculation of the Meissner masses in Sec. II.4. As the quark-gluon vertices in the calculation of the Meissner masses, the electroweak vertices have a nontrivial structure in Nambu-Gorkov space,

$$\Gamma_{ud/us,\pm}^{\mu}=\frac{e\,V_{ud/us}}{2\sqrt{2}\sin\theta_{W}}\left(\begin{array}[]{cc}\gamma^{\mu}(1-\gamma^{5})\,\tau_{ud/us,\pm}&0\\ 0&-\gamma^{\mu}(1+\gamma^{5})\,\tau_{ud/us,\mp}\end{array}\right)\,,$$ (128)

where $V_{ud}\simeq 0.97$ and $V_{us}\simeq 0.22$ are elements of the Cabibbo-Kobayashi-Maskawa (CKM) matrix, $\theta_{W}$ is the Weinberg angle, and

$$\tau_{ud,+} \equiv \left(\begin{array}[]{ccc}0&1&0\\ 0&0&0\\ 0&0&0\end{array}\right)\,,\qquad\tau_{ud,-}\equiv\left(\begin{array}[]{ccc}0&0&0\\ 1&0&0\\ 0&0&0\end{array}\right)\,,$$ (135)

$$\tau_{us,+} \equiv \left(\begin{array}[]{ccc}0&0&1\\ 0&0&0\\ 0&0&0\end{array}\right)\,,\qquad\tau_{us,-}\equiv\left(\begin{array}[]{ccc}0&0&0\\ 0&0&0\\ 1&0&0\end{array}\right)$$ (142)

are matrices in flavor space, with rows and columns labeled in the order $u,d,s$.

Finally, we have to specify the quark propagators. The structure in Nambu-Gorkov space from Eq. (20) carries over to the greater and lesser propagators,

$$S^{<,>}=\left(\begin{array}[]{cc}G_{+}^{<,>}&F_{-}^{<,>}\\[8.61108pt] F_{+}^{<,>}&G_{-}^{<,>}\end{array}\right)\,.$$ (143)

The normal propagators $G_{\pm}^{<,>}$ and anomalous propagators $F_{\pm}^{<,>}$ contain the color-flavor structure of the 2SC phase. From Eqs. (17) and (18) we read off the 2SC order parameter ${\cal M}=J_{3}I_{3}$, with $J_{A}$ and $I_{B}$ defined below Eq. (17). Therefore, ${\cal M}^{2}=J_{3}^{2}I_{3}^{2}$ with $J_{3}^{2}={\rm diag}(1,1,0)$ in color space and $I_{3}^{2}={\rm diag}(1,1,0)$ in flavor space. This shows that there are 4 gapped quasiparticles with identical gap and 5 ungapped quasiparticles. The corresponding projection operator for the gapped branches is simply ${\cal P}_{1}={\cal M}^{2}$. This structure shows that $ru$, $rd$, $gu$, and $gd$ quarks participate in pairing. For simplicity, we shall assume all quarks to be massless and the chemical potentials for $u$ and $d$ quarks to be identical,

$$\mu\equiv\mu_{u}=\mu_{d}\,.$$ (144)

Since we are interested in a deviation from chemical equilibrium, the chemical potential for strange quarks $\mu_{s}$ of course must be allowed to be different. Because of this distinction between unpaired strange quarks and unpaired up and down quarks, we introduce two projectors onto the unpaired sector, namely ${\cal P}_{2}=I_{3}^{2}(1-J_{3}^{2})$ (for $bu$ and $bd$ quarks) and ${\cal P}_{3}=1-I_{3}^{2}$ (for $rs$, $gs$, and $bs$ quarks). Obviously, ${\cal P}_{1}+{\cal P}_{2}+{\cal P}_{3}=1$. Then, we can write the normal and anomalous propagators as

$$G_{\pm}^{<,>}(K) = \sum_{r=1}^{3}{\cal P}_{r}G_{\pm,r}^{<,>}(K)\,\gamma^{0}\Lambda_{k}^{\mp}\,,$$ (145a)

$$F_{\pm}^{<,>}(K) = J_{3}I_{3}F_{\pm,1}^{<,>}(K)\,\gamma^{5}\Lambda_{k}^{\mp}\,,$$ (145b)

with the energy projectors $\Lambda_{k}^{\mp}$ from Eq. (10). We have omitted the anti-particle contribution (the $\pm$ signs correspond to charge conjugation). The anomalous propagators are proportional to the order parameter matrix ${\cal M}=J_{3}I_{3}$ and only the gapped branch, labeled by 1, appears. Written in this way, we have made all color, flavor, and Dirac structure explicit, and introduced greater and lesser propagators and their charge-conjugate counterparts for the three different branches $r=1,2,3$,

$$G_{\pm,r}^{>}(K) = -2\pi i\,\left\{B_{k,r}^{\pm}\,f(\epsilon_{k,r})\,\delta(k_{0}\pm\mu_{r}-\epsilon_{k,r})+B_{k,r}^{\mp}\,[1-f(\epsilon_{k,r})]\,\delta(k_{0}\pm\mu_{r}+\epsilon_{k,r})\right\}\,,$$ (146a)

$$G_{\pm,r}^{<}(K) = -2\pi i\,\left\{B_{k,r}^{\pm}\,[1-f(\epsilon_{k,r})]\,\delta(k_{0}\pm\mu_{r}-\epsilon_{k,r})+B_{k,r}^{\mp}\,f(\epsilon_{k,r})\,\delta(k_{0}\pm\mu_{r}+\epsilon_{k,r})\right\}\,,$$ (146b)

$$F_{\pm,r}^{>}(K) = 2\pi i\,\frac{\Delta}{2\epsilon_{k,r}}\Big\{f(\epsilon_{k,r})\,\delta(k_{0}\mp\mu_{r}-\epsilon_{k,r})-[1-f(\epsilon_{k,r})]\,\delta(k_{0}\mp\mu_{r}+\epsilon_{k,r})\Big\}\,,$$ (146c)

$$F_{\pm,r}^{<}(K) = 2\pi i\,\frac{\Delta}{2\epsilon_{k,r}}\Big\{[1-f(\epsilon_{k,r})]\,\delta(k_{0}\pm\mu_{r}-\epsilon_{k,r})-f(\epsilon_{k,r})\,\delta(k_{0}\pm\mu_{r}+\epsilon_{k,r})\Big\}\,,$$ (146d)

with $\mu_{1}=\mu_{2}=\mu$, $\mu_{3}=\mu_{s}$, the Fermi-Dirac distribution function $f(x)=1/(e^{x/T}+1)$, the dispersion relations

$$\epsilon_{k,1} = \sqrt{(k-\mu)^{2}+\Delta^{2}}\,,$$ (147a)

$$\epsilon_{k,2} = |k-\mu|\,,$$ (147b)

$$\epsilon_{k,3} = |k-\mu_{s}|\,,$$ (147c)

and the Bogoliubov coefficients

$$B_{k,r}^{\pm}\equiv\frac{\epsilon_{k,r}\pm(\mu_{r}-k)}{2\epsilon_{k,r}}\,.$$ (148)

These coefficients, and the existence of two terms in each of the propagators (146) reflects the fact that due to pairing the quasiparticles are mixtures of particles and holes, as shown in Fig. 3.

IV.3.2 $W$-boson polarization tensor

Having collected all vertices and propagators, we can now compute the $W$-boson polarization tensor. Inserting Eqs. (128) and (143) into Eq. (127) and performing the trace over Nambu-Gorkov space yields

$$\Pi_{\mu\nu}^{<,>}(Q) = -i\frac{e^{2}V_{us}^{2}}{8\sin^{2}\theta_{W}}\int_{P_{2}}\Big\{{\rm Tr}[\gamma^{\mu}(1-\gamma^{5})\,\tau_{us,+}\,G_{+}^{>,<}(P_{3})\,\gamma^{\nu}(1-\gamma^{5})\,\tau_{us,-}\,G_{+}^{<,>}(P_{2})]$$ (149)

$$+\;{\rm Tr}[\gamma^{\mu}(1+\gamma^{5})\,\tau_{us,-}\,G_{-}^{>,<}(P_{3})\,\gamma^{\nu}(1+\gamma^{5})\,\tau_{us,+}\,G_{-}^{<,>}(P_{2})]$$

$$-\;{\rm Tr}[\gamma^{\mu}(1-\gamma^{5})\,\tau_{us,+}\,F_{-}^{>,<}(P_{3})\,\gamma^{\nu}(1+\gamma^{5})\,\tau_{us,+}\,F_{+}^{<,>}(P_{2})]$$

$$-\;{\rm Tr}[\gamma^{\mu}(1+\gamma^{5})\,\tau_{us,-}\,F_{+}^{>,<}(P_{3})\,\gamma^{\nu}(1-\gamma^{5})\,\tau_{us,-}\,F_{-}^{<,>}(P_{2})]\Big\}\,.$$

For the color-flavor traces we need

$$\tau_{us,-}{\cal P}_{1} = \tau_{us,-}J_{3}^{2}\,,$$ (150a)

$$\tau_{us,-}{\cal P}_{2} = \tau_{us,-}(1-J_{3}^{2})\,,$$ (150b)

$$\tau_{us,+}{\cal P}_{3} = \tau_{us,+}\,,$$ (150c)

$$\tau_{us,+}{\cal P}_{1} = \tau_{us,+}{\cal P}_{2}=\tau_{us,-}{\cal P}_{3}=\tau_{us,\pm}I_{3}\tau_{us,\pm}I_{3}=0\,.$$ (150d)

In particular, this renders the contribution from the anomalous propagators zero. One can show that both terms in Eq. (149) with normal propagators yield the same contribution, and thus we can continue with twice the first term. Then, denoting

$${\cal T}^{\mu\nu}_{kq}\equiv{\rm Tr}[\gamma^{\mu}(1-\gamma^{5})\gamma^{0}\Lambda_{k}^{-}\gamma^{\nu}(1-\gamma^{5})\gamma^{0}\Lambda_{q}^{-}]\,,$$ (151)

we obtain

$$\Pi_{\mu\nu}^{<,>}(Q) = -i\frac{e^{2}V_{us}^{2}}{4\sin^{2}\theta_{W}}\int_{P_{2}}{\cal T}^{\mu\nu}_{p_{3}p_{2}}\left[2G_{+,3}^{>,<}(P_{3})G_{+,1}^{<,>}(P_{2})+G_{+,3}^{>,<}(P_{3})G_{+,2}^{<,>}(P_{2})\right]\,\,.$$ (152)

Here, the factor 2 in the first term is due to the color degeneracy and arises from the color trace over $J_{3}^{2}$. Inserting the propagators from Eqs. (146), using $1-f(\epsilon)=f(-\epsilon)$, and performing the integral over $u$ quark energies $p_{20}$ yields

$$\Pi_{\mu\nu}^{<,>}(Q) = i\frac{\pi e^{2}V_{us}^{2}}{2\sin^{2}\theta_{W}}\sum_{e_{2}e_{3}}\int_{\bm{p}_{2}}{\cal T}^{\mu\nu}_{p_{3}p_{2}}\Big\{2B_{p_{3},3}^{e_{3}}B_{p_{2},1}^{e_{2}}\,f(\pm e_{3}\epsilon_{p_{3},3})f(\mp e_{2}\epsilon_{p_{2},1})\,\delta(q_{0}-e_{3}\epsilon_{p_{3},3}+e_{2}\epsilon_{p_{2},1}+\delta\mu)$$ (153)

$$+\;B_{p_{3},3}^{e_{3}}B_{p_{2},2}^{e_{2}}\,f(\pm e_{3}\epsilon_{p_{3},3})f(\mp e_{2}\epsilon_{p_{2},2})\,\delta(q_{0}-e_{3}\epsilon_{p_{3},3}+e_{2}\epsilon_{p_{2},2}+\delta\mu)\Big\}\,,\hskip 28.45274pt$$

where the sums are over $e_{2},e_{3}=\pm 1$, where the upper (lower) sign corresponds to $\Pi^{<}$ ($\Pi^{>}$), and where we have abbreviated integration over three-momentum by

$$\int_{\bm{p}}\equiv\int\frac{d^{3}\bm{p}}{(2\pi)^{3}}\,.$$ (154)

IV.3.3 Collision integral

We now turn to evaluating the right-hand side of the kinetic equation (124), the ‘collision integral’. With the help of Eq. (125), we have

$${\rm Tr}[S^{>,<}(P_{1})\Sigma^{<,>}(P_{1})]=\frac{i}{M_{W}^{4}}\int_{P_{4}}{\rm Tr}[S^{>,<}(P_{1})\Gamma_{ud,-}^{\mu}S^{<,>}(P_{4})\Gamma_{ud,+}^{\nu}]\Pi^{>,<}_{\mu\nu}(Q)\,.$$ (155)

Taking the trace over Nambu Gorkov space yields

$${\rm Tr}[S^{>,<}(P_{1})\Gamma_{ud,-}^{\mu}S^{<,>}(P_{4})\Gamma_{ud,+}^{\nu}] = \frac{e^{2}V_{ud}^{2}}{8\sin^{2}\theta_{W}}\Big\{{\rm Tr}[G_{+}^{>,<}(P_{1})\gamma^{\mu}(1-\gamma^{5})\tau_{ud,-}G_{+}^{<,>}(P_{4})\gamma^{\nu}(1-\gamma^{5})\tau_{ud,+}]$$ (156)

$$+\;{\rm Tr}[G_{-}^{>,<}(P_{1})\gamma^{\mu}(1+\gamma^{5})\tau_{ud,+}G_{-}^{<,>}(P_{4})\gamma^{\nu}(1+\gamma^{5})\tau_{ud,-}]$$

$$-\;{\rm Tr}[F_{-}^{>,<}(P_{1})\gamma^{\mu}(1+\gamma^{5})\tau_{ud,+}F_{+}^{<,>}(P_{4})\gamma^{\nu}(1-\gamma^{5})\tau_{ud,+}]$$

$$-\;{\rm Tr}[F_{+}^{>,<}(P_{1})\gamma^{\mu}(1-\gamma^{5})\tau_{ud,-}F_{-}^{<,>}(P_{4})\gamma^{\nu}(1+\gamma^{5})\tau_{ud,-}]\Big\}\,.$$

For the traces over color and flavor space we need

$$\tau_{ud,\pm}{\cal P}_{1} = \tau_{ud,\pm}J_{3}^{2}\,,$$ (157a)

$$\tau_{ud,\pm}{\cal P}_{2} = \tau_{ud,\pm}(1-J_{3}^{2})\,,$$ (157b)

$$\tau_{ud,\pm}{\cal P}_{3} = 0\,,$$ (157c)

$$\tau_{ud,+}I_{3}\tau_{ud,+}I_{3} = I_{1}^{2}-1\,,$$ (157d)

$$\tau_{ud,-}I_{3}\tau_{ud,-}I_{3} = I_{2}^{2}-1\,.$$ (157e)

In contrast to the $W$-boson polarization tensor, now the traces containing anomalous propagators do not vanish. The four terms in Eq. (156) yield pairwise identical contributions, such that we may keep one of the terms containing normal propagators and one of the terms containing anomalous propagators and multiply the result by 2. We find

$${\rm Tr}[S^{>,<}(P_{1})\Sigma^{<,>}(P_{1})] = i\frac{e^{2}V_{ud}^{2}}{4M^{4}_{W}\sin^{2}\theta_{W}}\int_{P_{4}}\Big\{\left[2G_{+,1}^{>,<}(P_{1})G_{+,1}^{<,>}(P_{4})+G_{+,2}^{>,<}(P_{1})G_{+,2}^{<,>}(P_{4})\right]{\cal T}^{\mu\nu}_{p_{4}p_{1}}$$ (158)

$$+\;2F_{+,1}^{>,<}(P_{1})F_{-,1}^{<,>}(P_{4})\,{\cal U}^{\mu\nu}_{p_{4}p_{1}}\Big\}\,\Pi_{\mu\nu}^{>,<}(Q)\,,$$

where we have abbreviated the Dirac trace

$${\cal U}^{\mu\nu}_{kq}\equiv{\rm Tr}[\gamma^{\mu}(1-\gamma^{5})\gamma^{5}\Lambda_{k}^{+}\gamma^{\nu}(1+\gamma^{5})\gamma^{5}\Lambda_{q}^{-}]\,.$$ (159)

The different contributions in curly brackets in Eq. (158) describe the possible variants of the process $d\leftrightarrow W^{-}+u$ and can be understood as follows. The gapped $u$ and $d$ quarks, $r=1$, yield normal and anomalous contributions. Since there are two colors of gapped $u$ and $d$ quarks, these contributions come with a factor 2. The other contribution comes from the mode, $r=2$, the ungapped $u$ and $d$ quarks of a single color, which, of course, do not yield an anomalous contribution. Due to color conservation of the electroweak interaction, no mixing of the two modes is possible.

Eventually, we shall integrate both sides of the kinetic equation over the $d$-quark four-momentum $P_{1}=(p_{10},\bm{p}_{1})$ because we are interested in the total rate, not in the rate as a function of momentum. To this end, we need the result

$$\int_{-\infty}^{\infty}\frac{dp_{10}}{2\pi}{\rm Tr}[S^{>,<}(P_{1})\Sigma^{<,>}(P_{1})] = -i\frac{e^{2}V_{ud}^{2}}{4M_{W}^{4}\sin^{2}\theta_{W}}\sum_{e_{1}e_{4}}\int_{\bm{p}_{4}}\Bigg\{2\left[B_{p_{1},1}^{e_{1}}B_{p_{4},1}^{e_{4}}\,{\cal T}^{\mu\nu}_{p_{4}p_{1}}+\frac{e_{1}e_{4}\Delta^{2}}{4\epsilon_{p_{1},1}\epsilon_{p_{4},1}}{\cal U}^{\mu\nu}_{p_{4}p_{1}}\right]$$ (160)

$$\times f(\pm e_{1}\epsilon_{p_{1},1})f(\mp e_{4}\epsilon_{p_{4},1})\,\Pi_{\mu\nu}^{>,<}(e_{1}\epsilon_{p_{1},1}-e_{4}\epsilon_{p_{4},1},\bm{q})$$

$$+\,B_{p_{1},2}^{e_{1}}B_{p_{4},2}^{e_{4}}\,{\cal T}^{\mu\nu}_{p_{4}p_{1}}\,f(\pm e_{1}\epsilon_{p_{1},2})f(\mp e_{4}\epsilon_{p_{4},2})\,\Pi_{\mu\nu}^{>,<}(e_{1}\epsilon_{p_{1},2}-e_{4}\epsilon_{p_{4},2},\bm{q})\Bigg\}\,,$$

where we have performed the integration over the $u$ quark energy $p_{40}$ and the $d$ quark energy $p_{10}$, and where we have used the propagators (146). Now we insert the polarization tensors from Eq. (153) into Eq. (160). Evaluating the result requires us to calculate the contractions

$${\cal T}^{\mu\nu}_{p_{4}p_{1}}{\cal T}_{\mu\nu}^{p_{3}p_{2}} = 16(1-x_{34})(1-x_{12})\,,$$ (161a)

$${\cal U}^{\mu\nu}_{p_{4}p_{1}}{\cal T}_{\mu\nu}^{p_{3}p_{2}} = 8\left[(1-x_{34})(1-x_{12})+(1+x_{13})(1+x_{24})-(1+x_{14})(1+x_{23})\right]$$ (161b)

$$+8i\left[(\hat{\bm{p}}_{2}+\hat{\bm{p}}_{3})\cdot(\hat{\bm{p}}_{1}\times\hat{\bm{p}}_{4})+(\hat{\bm{p}}_{1}+\hat{\bm{p}}_{4})\cdot(\hat{\bm{p}}_{2}\times\hat{\bm{p}}_{3})\right]\,,\hskip 8.5359pt$$

where we have abbreviated

$$x_{mn}\equiv\hat{\bm{p}}_{m}\cdot\hat{\bm{p}}_{n}=\cos\theta_{mn}\,,$$ (162)

where $\theta_{mn}$ is the angle between the two three-momenta $\bm{p}_{m}$ and $\bm{p}_{n}$ with $m,n\in\{1,2,3,4\}$. Upon angular integration in the collision integral, the imaginary part in Eq. (161b) vanishes. Therefore, we drop this term from now on.

As mentioned above, the left-hand side of the kinetic equation (124) essentially becomes $\Gamma$, the net change of $d$ quarks per time and volume. In order to define $\Gamma$ properly and extract the correct numerical factor, we integrate the left-hand side over $P_{1}$ and project onto the $d$-flavor subspace. The latter is done by inserting the projector ${\cal Q}_{d}\equiv{\rm diag}(0,1,0)$. Moreover, to obtain the $d$ quark number density $n_{d}$ we also have to insert the matrix $\tau_{3}\equiv{\rm diag}(1,-1)$ in Nambu-Gorkov space (without pairing, the Nambu-Gorkov space would simply be a doubling of degrees of freedom, and the minus sign ensures that the charge-conjugate fermions yield the same contribution to the density as the original fermions). Then,

$$i\frac{\partial}{\partial t}\int_{P_{1}}{\rm Tr}[\gamma_{0}{\cal Q}_{d}\tau_{3}S^{<}(P_{1})]=-4\frac{\partial n_{d}}{\partial t}\equiv-4\Gamma\,,$$ (163)

where the factor 4 originates from the Nambu-Gorkov doubling and spin degeneracy.

Putting the results for the left-hand and right-hand sides of the kinetic equation together and dividing both sides by 4 yields

$$\Gamma=4\Gamma^{1131}+2\Gamma^{1231}+2\Gamma^{2132}+\Gamma^{2232}+4\widetilde{\Gamma}^{1131}+2\widetilde{\Gamma}^{1231}\,.$$ (164)

There are contributions from normal propagators,

$$\Gamma^{r_{1}r_{2}r_{3}r_{4}}\equiv 128\pi^{4}G_{F}^{2}V^{2}_{ud}V^{2}_{us}\sum_{e_{1}e_{2}e_{3}e_{4}}\int_{\bm{p}_{1}}\int_{\bm{p}_{2}}\int_{\bm{p}_{3}}\int_{\bm{p}_{4}}\,\delta(\bm{p}_{1}+\bm{p}_{2}-\bm{p}_{3}-\bm{p}_{4})\delta(e_{1}\epsilon_{1}+e_{2}\epsilon_{2}-e_{3}\epsilon_{3}-e_{4}\epsilon_{4}+\delta\mu)$$

$$\times\,(1-x_{12})(1-x_{34})\,B^{e_{1}}_{1}B^{e_{2}}_{2}B^{e_{3}}_{3}B^{e_{4}}_{4}\left[f(e_{1}\epsilon_{1})f(e_{2}\epsilon_{2})f(-e_{3}\epsilon_{3})f(-e_{4}\epsilon_{4})-f(-e_{1}\epsilon_{1})f(-e_{2}\epsilon_{2})f(e_{3}\epsilon_{3})f(e_{4}\epsilon_{4})\right]\,,\hskip 14.22636pt$$ (165)

and from anomalous propagators,

$$\widetilde{\Gamma}^{r_{1}r_{2}r_{3}r_{4}} \equiv 64\pi^{4}G_{F}^{2}V^{2}_{ud}V^{2}_{us}\sum_{e_{1}e_{2}e_{3}e_{4}}\int_{\bm{p}_{1}}\int_{\bm{p}_{2}}\int_{\bm{p}_{3}}\int_{\bm{p}_{4}}\,\delta(\bm{p}_{1}+\bm{p}_{2}-\bm{p}_{3}-\bm{p}_{4})\delta(e_{1}\epsilon_{1}+e_{2}\epsilon_{2}-e_{3}\epsilon_{3}-e_{4}\epsilon_{4}+\delta\mu)$$ (166)

$$\times\,\left[(1-x_{12})(1-x_{34})+(1+x_{13})(1+x_{24})-(1+x_{14})(1+x_{23})\right]\frac{e_{1}\Delta}{2\epsilon_{1}}B_{2}^{e_{2}}B_{3}^{e_{3}}\frac{e_{4}\Delta}{2\epsilon_{4}}$$

$$\times\,\left[f(e_{1}\epsilon_{1})f(e_{2}\epsilon_{2})f(-e_{3}\epsilon_{3})f(-e_{4}\epsilon_{4})-f(-e_{1}\epsilon_{1})f(-e_{2}\epsilon_{2})f(e_{3}\epsilon_{3})f(e_{4}\epsilon_{4})\right]\,.$$

We have expressed the factor from the electroweak vertices in terms of the Fermi coupling constant

$$G_{F}=\frac{\sqrt{2}e^{2}}{8M_{W}^{2}\sin^{2}\theta_{W}}=1.16637\cdot 10^{-11}\,{\rm MeV}^{-2}\,.$$ (167)

Also, we have introduced an additional integration over $\bm{p}_{3}$ and the corresponding $\delta$-function for momentum conservation, which was implicitly present in the previous expressions. Moreover, we have abbreviated the quasiparticle energies and Bogoliubov coefficients by

$$\epsilon_{i}\equiv\epsilon_{p_{i},r_{i}}\,,\qquad B_{i}^{e_{i}}\equiv B_{p_{i},r_{i}}^{e_{i}}\,.$$ (168)

Next we need to perform the angular integral. The 2SC phase is isotropic, i.e., none of the excitation energies and Bogoliubov coefficients depend on any angle, and thus we need to compute integrals of the form

$$K\equiv\int d\Omega_{1}\int d\Omega_{2}\int d\Omega_{3}\int d\Omega_{4}(1-x_{12})(1-x_{34})\delta(\bm{p}_{1}+\bm{p}_{2}-\bm{p}_{3}-\bm{p}_{4})\,.$$ (169)

We start with the $d\Omega_{2}$ integration and, to this end, abbreviate $\bm{P}\equiv\bm{p}_{3}+\bm{p}_{4}$ and factorize the $\delta$-function in terms of modulus and angles,

$$\delta[\bm{p}_{2}-(\bm{P}-\bm{p}_{1})]=\frac{1}{p_{2}^{2}}\delta(p_{2}-|\bm{P}-\bm{p}_{1}|)\delta(\Omega_{2}-\Omega)\,,$$ (170)

where $\Omega$ stands for the angles of the vector $\bm{P}-\bm{p}_{1}$. Thus, the effect of the angular $\delta$-function is to replace the direction $\hat{\bm{p}}_{2}$ with $(\bm{P}-\bm{p}_{1})/|\bm{P}-\bm{p}_{1}|$. Consequently, recalling the definition (162), we obtain

$$K=\frac{1}{p_{2}^{2}}\int d\Omega_{1}\int d\Omega_{3}\int d\Omega_{4}\left(1-\frac{\hat{\bm{p}}_{1}\cdot\bm{P}-p_{1}}{p_{2}}\right)(1-x_{34})\delta(p_{2}-|\bm{P}-\bm{p}_{1}|)\,.$$ (171)

Now, the $d\Omega_{1}$ integration can be performed by choosing the $z$-axis to be parallel to $\bm{P}$. We compute

$$\int d\Omega_{1}\left(1-\frac{\hat{\bm{p}}_{1}\cdot\bm{P}-p_{1}}{p_{2}}\right)\delta(p_{2}-|\bm{P}-\bm{p}_{1}|)=\frac{\pi}{p_{1}^{2}}\frac{1}{P}\Theta(P-p_{12})\Theta(P_{12}-P)(P_{12}^{2}-P^{2})\,,$$ (172)

where we have used $\Theta(a-|b|)=\Theta(a-b)\Theta(a+b)$ and introduced the abbreviations

$$p_{ij}\equiv|p_{i}-p_{j}|\,,\qquad P_{ij}\equiv p_{i}+p_{j}\,.$$ (173)

The remaining angular dependence is reduced to the angle between $\bm{p}_{3}$ and $\bm{p}_{4}$, which is present in $P=|\bm{p}_{3}+\bm{p}_{4}|$. Hence, one integration, say $d\Omega_{4}$, is trivial while we are left with the azimuthal integral from $d\Omega_{3}$,

$$K = \frac{8\pi^{3}}{p_{1}^{2}p_{2}^{2}}\int_{-1}^{1}dx\,\frac{1-x}{P}\Theta(P-p_{12})\Theta(P_{12}-P)(P_{12}^{2}-P^{2})$$ (174)

$$= \frac{4\pi^{3}}{p_{1}^{2}p_{2}^{2}p_{3}^{2}p_{4}^{2}}\int_{p_{34}}^{P_{34}}dP\,\Theta(P-p_{12})\Theta(P_{12}-P)(P_{12}^{2}-P^{2})(P_{34}^{2}-P^{2})\,.\hskip 28.45274pt$$

Since $P_{ij}>p_{ij}$, there are four orders of $p_{12},p_{34},P_{12},P_{34}$ that yield non-vanishing integrals,

$$p_{12}<p_{34}<P_{12}<P_{34}\,,$$ (175a)

$$p_{34}<p_{12}<P_{12}<P_{34}\,,$$ (175b)

$$p_{34}<p_{12}<P_{34}<P_{12}\,,$$ (175c)

$$p_{12}<p_{34}<P_{34}<P_{12}\,.$$ (175d)

The two other possibilities, $p_{12}<P_{12}<p_{34}<P_{34}$, and $p_{34}<P_{34}<p_{12}<P_{12}$ lead to vanishing integrals. Consequently,

$$K=\frac{4\pi^{3}}{p_{1}^{2}p_{2}^{2}p_{3}^{2}p_{4}^{2}}\,L(p_{12},P_{12},p_{34},P_{34})\,,$$ (176)

where we have defined

$$L(a,b,c,d) \equiv \Theta(c-a)\Theta(d-b)\Theta(b-c)\,J(c,b,b,d)+\Theta(a-c)\Theta(d-b)\,J(a,b,b,d)$$ (177)

$$+\,\Theta(a-c)\Theta(b-d)\Theta(d-a)\,J(a,d,b,d)+\Theta(c-a)\Theta(b-d)\,J(c,d,b,d)\,,$$

with

$$J(a,b,c,d) \equiv \int_{a}^{b}dP\,(c^{2}-P^{2})(d^{2}-P^{2})=c^{2}d^{2}(b-a)-\frac{1}{3}(c^{2}+d^{2})(b^{3}-a^{3})+\frac{1}{5}(b^{5}-a^{5})\,.$$ (178)

This is the result for the angular integration in Eq. (IV.3.3) and for the first term in Eq. (166). The two additional terms in Eq. (166) are easily obtained by a relabeling of momenta. Then, with the definitions

$$I(p_{1},p_{2},p_{3},p_{4}) \equiv L(p_{12},P_{12},p_{34},P_{34})\,,$$ (179a)

$$\tilde{I}(p_{1},p_{2},p_{3},p_{4}) \equiv L(p_{12},P_{12},p_{34},P_{34})+L(p_{24},P_{24},p_{13},P_{13})+L(p_{14},P_{14},p_{23},P_{23})\,,$$ (179b)

we arrive at

$$\Gamma^{r_{1}r_{2}r_{3}r_{4}} \equiv \frac{G_{F}^{2}V^{2}_{ud}V^{2}_{us}}{8\pi^{5}}\sum_{e_{1}e_{2}e_{3}e_{4}}\int_{p_{1}p_{2}p_{3}p_{4}}\,I(p_{1},p_{2},p_{3},p_{4})\,B^{e_{1}}_{1}B^{e_{2}}_{2}B^{e_{3}}_{3}B^{e_{4}}_{4}\delta(e_{1}\epsilon_{1}+e_{2}\epsilon_{2}-e_{3}\epsilon_{3}-e_{4}\epsilon_{4}+\delta\mu)$$ (180a)

$$\times\left[f(e_{1}\epsilon_{1})f(e_{2}\epsilon_{2})f(-e_{3}\epsilon_{3})f(-e_{4}\epsilon_{4})-f(-e_{1}\epsilon_{1})f(-e_{2}\epsilon_{2})f(e_{3}\epsilon_{3})f(e_{4}\epsilon_{4})\right]\,,$$

$$\widetilde{\Gamma}_{d}^{r_{1}r_{2}r_{3}r_{4}} \equiv \frac{G_{F}^{2}V^{2}_{ud}V^{2}_{us}}{16\pi^{5}}\sum_{e_{1}e_{2}e_{3}e_{4}}\int_{p_{1}p_{2}p_{3}p_{4}}\,\tilde{I}(p_{1},p_{2},p_{3},p_{4})\,\frac{e_{1}\Delta}{2\epsilon_{1}}B_{2}^{e_{2}}B_{3}^{e_{3}}\frac{e_{4}\Delta}{2\epsilon_{4}}\delta(e_{1}\epsilon_{1}+e_{2}\epsilon_{2}-e_{3}\epsilon_{3}-e_{4}\epsilon_{4}+\delta\mu)$$ (180b)

$$\times\left[f(e_{1}\epsilon_{1})f(e_{2}\epsilon_{2})f(-e_{3}\epsilon_{3})f(-e_{4}\epsilon_{4})-f(-e_{1}\epsilon_{1})f(-e_{2}\epsilon_{2})f(e_{3}\epsilon_{3})f(e_{4}\epsilon_{4})\right]\,,$$

where we have abbreviated

$$\int_{p_{1}p_{2}p_{3}p_{4}}\equiv\int_{0}^{\infty}dp_{1}\int_{0}^{\infty}dp_{2}\int_{0}^{\infty}dp_{3}\int_{0}^{\infty}dp_{4}\,.$$ (181)

Let us pause at this point in the calculation to explain the structure of the result we have obtained. We should be looking at Eqs. (164) and (180). In Eq. (164) we see that the rate receives contributions from 6 different terms. They correspond to 4 contributions where only normal propagators appear and 2 contributions which contain anomalous propagators. The 4 normal contributions contain all possible combinations of quasiparticle branches that can appear in the process $u+d\leftrightarrow u+s$. Recall that we labeled the branches by 1 (gapped $u$ and $d$ quarks), 2 (ungapped $u$ and $d$ quarks), and 3 (strange quarks, which are always ungapped). Therefore, for instance the rate $\Gamma^{1131}$ comes from the ingoing gapped $d$ and $u$ quarks (first pair of indices) and outgoing $s$ and $u$ quarks, of which the $u$ quark is gapped (second pair of indices). (The inverse process is also included, i.e., with $s$ and $u$ ingoing and $d$ and $u$ ingoing.) This contribution comes with a degeneracy factor 4 because at each of the electroweak vertices there are two possible colors, namely red and green, corresponding to the two colors of gapped $u$ and $d$ quarks. In contrast, the second term $\Gamma^{1231}$ only has degeneracy 2 because now the incoming $u$ quark is ungapped (labeled by 2), which fixes the color of the $u+W^{-}\to s$ vertex to be blue. The terms where both ingoing and outgoing quarks of one of the vertices are gapped (this can only happen at one vertex since the other vertex always has an ungapped $s$ quark attached to it) have counterparts where both gapped quarks are represented by anomalous propagators, $\widetilde{\Gamma}^{1131}$ and $\widetilde{\Gamma}^{1231}$ in Eq. (164). We show all 6 terms in diagrammatic form in Fig. 8.

Figure 8: Subprocesses contributing to the rate of the process $u+d\to u+s$ in the 2SC phase. Gapped branches are indicated by a $\Delta$ at the quark lines, with $s$ quarks being always ungapped. The two last diagrams contain anomalous propagators, indicated by the condensate insertion.

Notice that one diagram (with degeneracy 1) only contains ungapped quarks. At low temperatures, the total rate is well approximated by this diagram because a gap in at least one of the dispersions of the participating quasiparticles induces exponential suppression at temperatures much smaller than the gap. This exponential suppression, however, can be alleviated by a sufficiently large $\delta\mu$. For a detailed discussion see Ref. Alford:2006gy.

The calculation of the rate is further complicated by the fact that each term in Eq. (164) is in turn composed of various subprocesses, as Eq. (180) shows. Due to the sum over $e_{1},e_{2},e_{3},e_{4}=\pm$ there are in fact 16 such subprocesses (not all of them actually giving a nonzero contribution). What is the meaning of these subprocesses? Naively, had we written down an ‘ordinary’ collision integral for the process $u+d\to u+s$ we would have expected a structure of the Fermi distribution functions of the form $f_{d}f_{u}(1-f_{s})(1-f_{u})$ because this corresponds to the probability of finding the ingoing states occupied and the outgoing states empty. However, with $f(-x)=1-f(x)$, Eq. (180) shows that all possible combinations of $f$ and $1-f$ formally appear. Therefore, the quasiparticles, which actually are momentum-dependent mixtures of particles and holes, can appear on either side of the reaction. This is a typical property of fermionic superconductors and superfluids, where particles can be created from or deposited into the condensate since particle number is spontaneously broken. The same effect is thus also present in superfluid ³He Vollhardt1990. For a more detailed study of the contribution of the various subprocesses, based on a numerical evaluation, see Ref. Alford:2006gy.

IV.3.4 Unpaired limit

It is instructive to take the limit of unpaired quark matter from our complicated result. To this end, it is useful to slightly rewrite the Bogoliubov coefficients and the dispersion relation. We define

$$\tilde{B}_{i}^{e_{i}}\equiv\frac{1}{2}\left(1+e_{i}\frac{p_{i}-\mu_{i}}{\tilde{\epsilon}_{i}}\right)\,,\qquad\tilde{\epsilon}_{i}\equiv{\rm sgn}(p_{i}-\mu_{i})\epsilon_{i}\,.$$ (182)

Then we note that for any function $F$ we have

$$\sum_{e_{i}}\int_{0}^{\infty}dp_{i}\,B_{i}^{e_{i}}F(e_{i}\epsilon_{i})=\sum_{e_{i}}\int_{0}^{\infty}dp_{i}\,\tilde{B}_{i}^{e_{i}}F(-e_{i}\tilde{\epsilon}_{i})\,.$$ (183)

Consequently, in Eqs. (180) we can simply replace $B_{i}^{e_{i}}$ by $\tilde{B}_{i}^{e_{i}}$ and $e_{i}\epsilon_{i}$ by $-e_{i}\tilde{\epsilon}_{i}$. The advantage of the new formulation is that now $\tilde{B}_{i}^{+}=1$ and $\tilde{B}_{i}^{-}=0$ if $\Delta=0$. More physically speaking, as we let $\Delta\to 0$, the dispersions $\tilde{\epsilon}_{i}$ and $-\tilde{\epsilon}_{i}$ reduce to pure particle and hole branches, respectively, while $\epsilon_{i}$ becomes a hole below the Fermi surface and a particle above the Fermi surface and vice versa for $-\epsilon_{i}$. Therefore, after the replacement with the new Bogoliubov coefficients and excitation energies, the rate of the unpaired phase is given by the term $e_{1}=e_{2}=e_{3}=e_{4}=+1$ in Eq. (180a) [the anomalous contribution (180b) obviously vanishes completely for $\Delta=0$]. We therefore obtain

$$\Gamma_{0}=\frac{N_{c}^{2}G_{F}^{2}V^{2}_{ud}V^{2}_{us}}{8\pi^{5}}\int_{p_{1}p_{2}p_{3}p_{4}}\,I(p_{1},p_{2},p_{3},p_{4})\,\delta(p_{1}+p_{2}-p_{3}-p_{4})$$

$$\times\,\Big\{[1-f(p_{1}-\mu)][1-f(p_{2}-\mu)]f(p_{3}-\mu_{s})f(p_{4}-\mu)-f(p_{1}-\mu)f(p_{2}-\mu)[1-f(p_{3}-\mu_{s})][1-f(p_{4}-\mu)]\Big\}\,,\hskip 19.91684pt$$ (184)

where we have reinstated the degeneracy factor $N_{c}^{2}$ – one factor of $N_{c}$ for each of the electroweak vertices – to obtain the complete result for three-flavor unpaired quark matter.

We can now make use of the $\delta$-function to perform the $p_{4}$ integral, which yields a factor $\Theta(p_{1}+p_{2}-p_{3})$ in the integrand. From now on, we shall assume $\delta\mu>0$. This means that the chemical potential for strange quarks is larger than in equilibrium and thus the system will react by converting $s$ quarks into $d$ quarks. Therefore, the first term in the curly brackets will give a larger contribution than the second term. It is thus convenient to use the explicit form of the Fermi distribution functions to rewrite the second term,

$$\Gamma_{0} = \frac{2N_{c}^{2}G_{F}^{2}V^{2}_{ud}V^{2}_{us}}{15\pi^{5}}\left(1-e^{-\delta\mu/T}\right)\int_{0}^{\infty}dp_{1}\int_{0}^{\infty}dp_{2}\int_{0}^{\infty}dp_{3}\,H(p_{1},p_{2},p_{3})$$ (185)

$$\times\,[1-f(p_{1}-\mu)][1-f(p_{2}-\mu)]f(p_{3}-\mu_{s})f(p_{1}+p_{2}-p_{3}-\mu)\,,$$

where we have also explicitly evaluated the function $I(p_{1},p_{2},p_{3},p_{4})$, which, using Eqs. (177), (178), and (179a), yields

$$H(p_{1},p_{2},p_{3}) \equiv [\Theta(p_{2}-p_{3})-\Theta(p_{1}-p_{3})]^{2}[\Theta(p_{1}-p_{2})h_{1}+\Theta(p_{2}-p_{1})h_{2}]$$ (186)

$$+[\Theta(p_{1}+p_{2}-p_{3})-\Theta(p_{2}-p_{3})-\Theta(p_{1}-p_{3})]\Theta(2p_{3}-p_{1}-p_{2})h_{3}$$

$$+\{1-[\Theta(p_{2}-p_{3})-\Theta(p_{1}-p_{3})]^{2}\}\Theta(p_{1}+p_{2}-2p_{3})h_{4}\,,$$

with

$$h_{1} \equiv p_{2}^{3}(10p_{1}^{2}+5p_{1}p_{2}+p_{2}^{2})\,,$$ (187a)

$$h_{2} \equiv p_{1}^{3}(10p_{2}^{2}+5p_{1}p_{2}+p_{1}^{2})\,,$$ (187b)

$$h_{3} \equiv (p_{1}+p_{2}-p_{3})^{3}[(p_{1}+p_{2})^{2}+3(p_{1}+p_{2})p_{3}+6p_{3}^{2}]\,,$$ (187c)

$$h_{4} \equiv p_{3}^{3}[10(p_{1}+p_{2})^{2}-15(p_{1}+p_{2})p_{3}+6p_{3}^{2}]\,.$$ (187d)

The integral in Eq. (185) is still too complicated to be evaluated analytically for arbitrary $\mu$, $T$, and $\delta\mu$. We shall therefore discuss two limit cases where analytical results can be obtained. Firstly, let us discuss the case $T=0$. In this case, the Fermi functions in Eq. (185) become step functions, $\lim_{T\to 0}f(x)=\Theta(-x)$. As a consequence, all step functions in Eq. (186) are 1 or 0 in the entire integration domain except for $\Theta(p_{1}+p_{2}-p_{3})$. Due to this step function, one has to distinguish the cases $\mu_{s}<2\mu$ and $\mu_{s}>2\mu$. We are mostly interested in small deviations from chemical equilibrium. Therefore, we omit the case $\mu_{s}>2\mu$ (the result can be found in Ref. PhysRevD.47.325). For the other case and setting $N_{c}=3$ we find

$$\Gamma_{0}(T=0,\delta\mu<\mu) = \frac{6}{5\pi^{5}}G_{F}^{2}V^{2}_{ud}V^{2}_{us}\int_{\mu}^{\mu_{s}}dp_{3}\int_{\mu}^{p_{3}}dp_{2}\int_{\mu}^{p_{3}-p_{2}+\mu}dp_{1}\,h_{3}$$ (188)

$$= \frac{16}{5\pi^{5}}G_{F}^{2}V_{us}^{2}V_{ud}^{2}\left(\mu^{5}\delta\mu^{3}+\frac{5}{16}\mu^{4}\delta\mu^{4}-\frac{3}{16}\mu^{3}\delta\mu^{5}+\frac{1}{32}\mu^{2}\delta\mu^{6}+\frac{5}{112}\mu\delta\mu^{7}-\frac{15}{896}\delta\mu^{8}\right)\,.$$

Secondly, we compute the rate for $\delta\mu\ll T\ll\mu$ (‘subthermal regime’). In this case we observe that the prefactor of the integral in Eq. (185) is linear in $\delta\mu/T$ to lowest order. Neglecting all higher order corrections, we can thus set $\delta\mu=0$ in the integral. We then introduce the dimensionless integration variables

$$x=\frac{p_{1}-\mu}{T}\,,\qquad y=\frac{p_{2}-\mu}{T}\,,\qquad z=\frac{p_{3}-\mu_{s}}{T}\,,$$ (189)

and approximate the lower boundaries of the $x$, $y$ and $z$ integrals by $-\infty$ because of $\mu,\mu_{s}\gg T$. Moreover, we may expand the functions $h_{i}$ for large $\mu/T$, which yields to lowest order $h_{1}\simeq h_{2}\simeq h_{3}\simeq h_{4}\simeq 16\mu^{5}$. Then, all step functions in Eq. (186) combine to 1 such that we arrive at

$$\Gamma_{0}(\delta\mu\ll T\ll\mu)=\frac{64G_{F}^{2}V_{ud}^{2}V_{us}^{2}}{5\pi^{3}}T^{2}\mu^{5}\delta\mu\,,$$ (190)

where we have used

$$\int_{-\infty}^{\infty}dx\int_{-\infty}^{\infty}dy\int_{-\infty}^{\infty}dz\,[1-f(x)][1-f(y)]f(z)f(x+y-z)=\frac{2\pi^{2}}{3}\,.$$ (191)

Here, in a slight abuse of notation, the Fermi-Dirac function is $f(x)=1/(e^{x}+1)$, not $f(x)=1/(e^{x/T}+1)$. The result (190) was first computed in Ref. Heiselberg:1992bd.

IV.4 Comparing different color superconductors

To summarize, we have computed the net $d$ quark production rate of the processes $u+d\leftrightarrow s+d$. We have derived the result for 2SC quark matter in the form of collision integrals which have to be evaluated numerically, and we have derived analytical limits for the case of unpaired quark matter. It remains to insert the result of the rate into the expression (123) for the frequency-dependent bulk viscosity. In this expression we already have assumed the rate to be linear in the deviation from equilibrium $\delta\mu$. Therefore, for unpaired quark matter, consistency requires us to work with the approximation (190). As this analytical result shows explicitly, the reaction rate increases monotonically (and quadratically) with temperature. Therefore, as already anticipated below Eq. (123), the bulk viscosity is expected to have a maximum for a certain temperature at given external frequency $\omega$. This behavior is borne out in Fig. 9, where the (red and green) curves for unpaired and 2SC quark matter are taken from Ref. Alford:2006gy. For essentially the entire temperature range shown in this figure, the rate of $u+d\leftrightarrow s+d$ in the 2SC phase is well approximated by 1/9 of the value of the unpaired phase, due to the exponential suppression of all sub-processes that involve one or more gapped excitation branches (the green 2SC curve in Fig. 9 assumes a critical temperature for the 2SC phase of $T_{c}=30\,{\rm MeV}$). While the electroweak rate of the 2SC phase is smaller than that of the unpaired phase, this is not true for the bulk viscosity. At a given temperature, it is not the largest reaction rate that gives the largest viscosity but the rate that is closest to the external frequency. For the slower rate of the 2SC phase a larger temperature is needed to be in resonance with the external oscillation. That is why the maximum of the bulk viscosity is shifted to a higher temperature due to Cooper pairing.

Figure 9: Bulk viscosity $\zeta$ as a function of temperature for various color superconductors and unpaired quark matter. For all curves the external oscillation frequency is $\omega/(2\pi)=1\,{\rm ms}^{-1}$, and for all curves except for the inset the quark chemical potential is $\mu=400\,{\rm MeV}$. The calculation presented in detail in the text gives the (red and green) curves for unpaired and 2SC results from the process $u+d\leftrightarrow u+s$ Alford:2006gy. The results for CFL are shown for thermal kaons Alford:2007rw, $\delta m\equiv m_{K^{0}}-\mu_{K^{0}}>0$, and in the presence of kaon condensation Alford:2008pb, $\delta m<0$. In CFL there is also a contribution purely from the superfluid phonon $\phi$ Manuel:2007pz. The inset shows the 2SC result from direct Urca processes Alford:2025tbp for a baryon density of 7 times (nuclear) saturation density, employing an NJL model with different ratios for vector and scalar couplings, $G_{V}/G_{S}=0.6,0.8,1.0,1.2$ from top to bottom.

The results of Fig. 9 require an input for the thermodynamic functions $B$ and $C$, defined in Eqs. (117). For the unpaired and 2SC curves from the process $u+d\leftrightarrow u+s$ shown here, the relations for non-interacting, zero-temperature quark matter have been used (augmented by the superconducting gap in the case of 2SC). As discussed above, see Eq. (121), it is crucial to include a nonzero strange quark mass in the equilibrium quantities $B$ and $C$ because otherwise compression and expansion cannot induce an out-of-equilibrium situation (in the 2SC phase, the gap alone would suffice for this purpose because strange quarks are ungapped). Including the effect of interactions into the thermodynamic functions yields corrections to the bulk viscosity. For unpaired quark matter, this question was addressed within perturbative QCD in Ref. CruzRojas:2024etx, where also a comparison to strong-coupling results from the gauge-gravity duality was discussed. In the inset of Fig. 9 we show 2SC results from direct Urca processes (112), where interactions computed in an NJL model are taken into account to evaluate the relevant susceptibilities. The curves shown in the inset, taken from Ref. Alford:2025tbp, are valid for a fixed density that does not correspond exactly to the chemical potential of the other curves, and four different curves are shown for different values of the coupling strength in the model. Therefore, a direct quantitative comparison to the other curves has to be taken with care, the main point is rather as follows: If two reactions contribute to the same out-of-equilibrium chemical potential $\delta\mu$, the faster one dominates the bulk viscosity, even if the bulk viscosity was larger in the absence of the fast process. Here, however, Urca processes and non-leptonic processes re-equilibrate different $\delta\mu$’s. Therefore, the contribution from the Urca process, whose rate is slower than that of the non-leptonic one, can enhance the bulk viscosity at large temperatures and the inset suggests that this contribution is indeed relevant for instance for neutron star mergers. In unpaired quark matter, both classes of reactions were included within a single calculation in Ref. Sad:2007afd.

The other main point of Fig. 9 is the comparison to the CFL bulk viscosity. In that case, the contribution of the quarks is negligible because all quarks acquire a gap in their dispersion. The main contribution comes from strangeness changing processes involving kaons. As discussed briefly in Sec. III, it is expected that the CFL kaons condense. The figure shows the results for both cases, i.e., uncondensed thermal Alford:2007rw and condensed kaons Alford:2008pb. In both cases, the bulk viscosity arises from the interaction of the kaons with the Goldstone mode from superfluidity. The figure also shows the contribution from processes involving the Goldstone mode alone Manuel:2007pz. We see that the bulk viscosity of CFL is, for most relevant temperatures, negligibly small compared to that of 2SC and unpaired quark matter. Since the electroweak rate from kaons is smaller than the one from fermions, the bulk viscosity peaks at much larger temperatures. We see that a peak starts to develop at $T>10\,{\rm MeV}$, a regime particularly interesting for neutron star mergers. However, we need to keep in mind that in this regime we are approaching the critical temperature of kaon condensation Alford:2007qa, and, in fact, the critical temperature of CFL itself.

V Open questions and future directions

Color superconductivity is well understood at ultra-high quark densities, where the QCD coupling is weak. Weak-coupling methods fail at moderate densities, and we have no first-principle understanding of color superconductivity in the regime relevant for neutron stars. Even if we were to ignore Cooper pairing or a strong coupling version thereof, the sign problem of lattice QCD prevents us from calculating by brute force the properties of cold and dense matter. Progress to overcome or circumvent the sign problem has been made from various different angles Aarts:2008rr; Cristoforetti:2014gsa; Langelage:2014vpa; Alexandru:2015sua; Marchis:2017oqi; Ratti:2022qgf; Basar:2023bwd, but further breakthrough ideas are needed to compute the QCD phase diagram at finite density on the lattice.

Besides the strong-coupling nature of the problem, color superconductivity at moderately high densities is difficult due to the effect of the strange quark mass, which disrupts the particularly symmetric paring pattern of CFL. A lot of progress has been made with the help of effective theories and phenomenological models, which give us an idea about the possible, less symmetric, candidate phases. It is still an open question, however, whether the CFL phase persists down to densities where the transition to nuclear matter occurs or whether there are non-CFL color superconductors in between. These phases are likely to contain unpaired fermionic modes, or at the very least fermionic modes with very small gaps, and possibly break rotational and/or translational invariance. One possible route to better understand the breakdown density of CFL is to build on the recent progress in dense perturbative QCD with unpaired quark matter Kurkela:2009gj; Gorda:2021znl; Fernandez:2021jfr; Gorda:2023mkk. Using these results and ideally combining them with the gap equation, progress could be made towards a more reliable (weak-coupling) calculation of the ratio $m_{s}^{2}/(\mu\Delta)$, which determines the mismatch in Fermi momenta that seeks to break CFL.

The main ‘laboratory’ to test the predictions and conjectured phases of color superconductivity are neutron stars. (It is conceivable that heavy-ion collisions may shed some light on color superconductivity in future experiments as well Kitazawa:2001ft; Schmitt:2016pre; Nishimura:2024kvz.) Neutron stars are of course complicated objects, and many of their observable properties depend not only on the core but also on the outer low-density layers. It is therefore very difficult to find a ‘smoking gun’ signature for the existence of deconfined quark matter (see Ref. Annala:2019puf for an indication based on the equation of state) or color superconductivity. Nevertheless, more, and more precise, astrophysical data has become available and will become available in the near future. Most notably, new insights can be expected from the detection of gravitational waves from neutron star mergers LIGOScientific:2018hze and possibly from isolated neutron stars 2024arXiv240302066J, with updated current Acernese_2015; Aasi_2015 and new Abbott_2017; Maggiore_2020 detectors, besides improved data to ‘measure’ neutron star masses and radii, which also has seen impressive progress in recent years NANOGrav:2019jur; Riley:2019yda; Miller:2019cac; Fonseca:2021wxt; Miller:2021qha; Riley:2021pdl; Romani:2022jhd; Choudhury:2024xbk.

It is an interesting question whether color superconductivity plays a detectable role in a neutron star merger. Estimates of critical temperatures of most color superconductors indicate that Cooper pairing might survive even in the hot environment of a merger process. Therefore, it is worth investigating, for instance, if superfluid properties of CFL need to be taken into account in the hydrodynamic simulation, potentially within a two-fluid formalism that distinguishes between superfluid and normal components. Another aspect is the dynamical transition from hadronic to quark matter in a merger process, for instance if two neutron stars without quark core merge and the phase transition is triggered by the creation of hot and dense spatial regions. Even without color superconductivity this transition is treated non-dynamically in current simulations Most:2019onn; Tootle:2022pvd; Guo:2023som. While adding color superconductivity to the equation of state is straightforward, further studies are needed to understand the effect of color superconductivity on the dynamical phase conversion under merger conditions. One may also ask whether color superconductivity can be linked to continuous gravitational waves from isolated neutron stars. Besides originating from potentially unstable $r$-modes, they can also be generated by ‘mountains’, non-axisymmetric deformations of the star sustained by the internal structure. While the physics of the crust is the main candidate in the literature for the source of mountains, color superconductivity in the core can play a role if it is in a crystalline phase Knippel:2009st or if it generates a lattice of magnetic flux tubes Glampedakis:2012qp; Haber:2017oqb.

Astrophysical constraints from mass and radius measurements and the tidal deformability extracted from neutron star mergers can be used for a model-independent constraint for the color-superconducting gap, as pointed out recently Kurkela:2024xfh; Kurkela:2025acm. The idea is to compute an upper bound on the condensation energy (72) from the astrophysical constraints, thermodynamic consistency, and the low-density nuclear matter equation of state Hebeler:2013nza. This approach can be further improved in an obvious way by new, more constraining, data, and by further improvements in our understanding of the low-density and ultra-high-density limits of the equation of state. The approach is agnostic to the mechanism that lowers the free energy of unpaired quark matter. One may ask whether other exotic non-perturbative phenomena compete with or coexist with color superconductivity at moderate densities. One candidate is the quarkyonic phase, the large-$N_{c}$ version of nuclear matter McLerran:2007qj, modeled with the help of a layered structure of the Fermi sphere McLerran:2018hbz. It is unknown whether a remnant or variant of this phase survives at $N_{c}=3$. Recently, Cooper pairing of quarks in the quarkyonic phase has been considered within a phenomenological approach Gartlein:2025zhd, and it remains to be seen whether this phase can be confirmed in a more microscopic treatment. The quarkyonic phase is obtained fully dynamically within the gauge-gravity duality Kovensky:2020xif. A natural question for future research is how realistic the holographic description of color superconductivity can be made. Two main lines of ‘holographic color superconductivity’ exist, either focusing on the spontaneous breaking of color symmetry Chen:2009kx; Faedo:2018fjw; Preau:2025ubr or on simply mimicking the Cooper pair condensate by a scalar field Basu:2011yg; BitaghsirFadafan:2018iqr; CruzRojas:2025fzs. In either case, the theories under consideration are different from QCD and further studies towards a more realistic strong-coupling picture of color superconductivity are needed.

Another unsolved fundamental question is the one about the quark-hadron continuity, a particular aspect of the QCD phase structure at high baryon densities. This question is often discussed without taking into account pairing, especially if the focus is on the equation of state and its astrophysical implications Baym:2017whm; Baym:2019iky. In that case, since neither order parameter for deconfinement and chiral transitions corresponds to an exact symmetry, a continuous transition from nuclear matter to quark matter is allowed. This is more complicated in the presence of pairing, where the exact symmetry of baryon number conservation can be broken spontaneously, resulting in a superfluid. CFL is a superfluid, and although the 2SC in its usual definition is not, additional single flavor pairing can render this phase a superfluid too Fujimoto:2019sxg. It was pointed out in Refs. Schafer:1998ef; Alford:1999pa that a continuity between (hyper)nuclear matter and CFL is possible, and the role of the axial anomaly for the presence of a crossover and a high-density critical point was discussed in Refs. Hatsuda:2006ps; Yamamoto:2007ah; Schmitt:2010pf. The question of the quark-hadron continuity is not only of theoretical interest for the QCD phase diagram, but may be realized at a quark-hadron interface inside neutron stars. A particularly interesting aspect of the interface between CFL and superfluid nuclear matter, yet to be fully understood, is the possible existence of continuous vortices Alford:2018mqj; Chatterjee:2018nxe, also considered in two-flavor quark matter Fujimoto:2021bes. Vortices are also, on the other hand, argued to provide a qualitative difference between quark and hadronic matter despite the identical global symmetries Cherman:2018jir; Cherman:2020hbe.

Astrophysical data (and possibly data from heavy-ion collisions) are expected to improve our understanding of color superconductivity, but the intricacies of different color-superconducting phases and the transition between them will always be difficult to extract from these highly complex systems. And, of course, we do not know for certain that the conditions in these systems are extreme enough to even exhibit color superconductivity. Therefore, one can ask if there is a more controlled (tabletop) experiment that can at least mimic color-superconducting matter. Experiments with ultracold fermionic atoms are a very good candidate. In the simplest situation of two fermionic species they show a BEC-BCS crossover 2008NCimR..31..247K. This crossover is conjectured in QCD at finite isospin chemical potential (and vanishing baryon chemical potential) Son:2000xc and in $N_{c}=2$ QCD Iida:2019rah; Boz:2019enj. However, in real-world dense QCD the situation is more complicated since we know that at low densities bound states are made of three quarks, not two. It is nevertheless desirable to reproduce certain aspects of color-superconducting quark matter in cold atomic systems, for instance pairing with mismatched Fermi surfaces Zwierlein:2006zz; Partridge:2006zkb. Also, atomic systems are routinely being prepared with more that two species and with ‘artificial’ (non-abelian) gauge fields PhysRevLett.98.160405; PhysRevLett.101.203202; PhysRevLett.102.130401; PhysRevLett.109.240401; PhysRevA.97.023632; Wang_2020; 2025arXiv250204714M and it would be very interesting to create, for instance, an atomic analogue of color-flavor locking.