# Topological susceptibility and axion potential in two-flavor superconductive
quark matter

Fabrizio Murgana<EMAIL_ADDRESS>Department of Physics and
Astronomy ”Ettore Majorana”, University of Catania, Via S. Sofia 64, I-95123
Catania, Italy  INFN-Sezione di Catania, Via S. Sofia 64, I-95123 Catania,
Italy  David E. Alvarez Castillo<EMAIL_ADDRESS>Institute of Nuclear
Physics, Polish Academy of Sciences, Radzikowskiego 152, 31-342 Cracow, Poland
Incubator of Scientific Excellence - Centre for Simulations of Superdense
Fluids, University of Wroclaw, plac Maksa Borna 9, PL-50204 Wroclaw, Poland
Facultad de Ciencias Físico Matemáticas, Universidad Autónoma de Nuevo León,
Av. Universidad S/N, C.U., 66455 San Nicolás de los Garza, N.L., Mexico  Ana
G. Grunfeld<EMAIL_ADDRESS>CONICET, Godoy Cruz 2290, C1425FQB
Ciudad Autónoma de Buenos Aires, Argentina  Departamento de Física, Comisión
Nacional de Energía Atómica, Avenida Libertador 8250, C1429 BNP, Ciudad
Autónoma de Buenos Aires, Argentina  Marco Ruggieri
<EMAIL_ADDRESS>Department of Physics and Astronomy ”Ettore
Majorana”, University of Catania, Via S. Sofia 64, I-95123 Catania, Italy
INFN-Sezione di Catania, Via S. Sofia 64, I-95123 Catania, Italy

###### Abstract

We study the potential of the axion, $a$, of Quantum Chromodynamics, in the
two-flavor color superconducting phase of cold and dense quark matter. We
adopt a Nambu-Jona-Lasinio-like model. Our interaction contains two terms, one
preserving and one breaking the $U(1)_{A}$ symmetry: the latter is responsible
of the coupling of axions to quarks. We introduce two quark condensates,
$h_{L}$ and $h_{R}$, describing condensation for left-handed and right-handed
quarks respectively; we then study the loci of the minima of the thermodynamic
potential, $\Omega$, in the $(h_{L},h_{R})$ plane, noticing how the instanton-
induced interaction favors condensation in the scalar channel when the
$\theta-$angle, $\theta=a/f_{a}$, vanishes. Increasing $\theta$ we find a
phase transition where the scalar condensate rotates into a pseudo-scalar one.
We present an analytical result for the topological susceptibility, $\chi$, in
the superconductive phase, which stands both at zero and at finite
temperature. Finally, we compute the axion mass and its self-coupling. In
particular, the axion mass $m_{a}$ is related to the full topological
susceptibility via $\chi=m_{a}^{2}f_{a}^{2}$, hence our result for $\chi$
gives an analytical result for $m_{a}$ in the superconductive phase of high-
density Quantum Chromodynamics.

QCD axion, QCD phase diagram, chiral symmetry restoration, compact stellar
objects

###### pacs:

12.38.Aw,12.38.Mh

##  I Introduction

Axions, originally proposed by Peccei and Quinn in 1977 as a possible solution
to the strong CP problem in Quantum Chromodynamics (QCD) [1, 2, 3, 4, 5, 6, 7,
8], have since then become prime candidates for dark matter [3, 9, 10, 11, 12,
13, 14, 15]. The quest to understand the nature of dark matter, which
constitutes a significant portion of the universe’s mass, has led to extensive
theoretical and experimental efforts, with axions emerging as particularly
compelling candidates due to their unique properties. Theoretical models and
experiments constrain the axion mass to be very light, on the order of
$10^{-6}$ to $10^{-3}$ eV [16]. This range arises from considerations such as
the possible violation of the CP symmetry of QCD and the consequent electric
dipole moment for the neutron [17, 18, 19, 20, 21, 22, 23, 24], astrophysical
and cosmological observations, and experimental searches for axions using
techniques like cavity haloscopes and axion helioscopes. In addition to their
potential role as dark matter candidates, axions, behaving as scalar fields,
may also serve as the constituents of boson stars [25] and arrange into axion
stars [26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40] and form
Bose-Einstein condensates [41, 42].

Furthermore, not only axions but axion-like particles (ALP) have been proposed
to play the role of dark matter in order to account for missing matter at
cosmological scales. It has been speculated that axions may exist in compact
star interiors, and can be buried deep in their cores, giving rise to
interaction with quark matter, as we investigate in this study. In the
assumption that axions may be able to escape compact star interiors, they
would be prone to interact with magnetic fields via the so called Primakov
effect, i.e., the axion resonantly converting into a radio photon [43, 44],
which would provide a feasible way to detect them. The emission of axions
could also cool compact stars, deviating their thermal evolution from standar
scenarios [45, 46, 47, 48].

Macroscopic properties of compact stars bearing axions have been presented in
[49, 50], where particular modeling has resulted in unstable compact stars
that would collapse due to radial oscillations if vector repulsive
interactions in quark matter are not taken into account. Furthermore, axions
have been investigated within astrophysical contexts, particularly in relation
to supernova explosions and the formation of protoneutron stars [51, 52].
Previous works on the coupling of the QCD axion to quarks can be found in [63,
53, 55, 56, 54, 57, 58, 59, 49, 60, 61, 62], where in particular the effect of
the QCD chiral phase transition on the low-energy properties of the axion
itself has been investigated, both at finite temperature and at finite quark
chemical potential, $\mu$. The axions enter the model similarly to the
$\theta-$angle of QCD; in fact, formally one can pass from QCD at finite
$\theta$ to QCD with a finite axion background by identifying
$\theta=a/f_{a}$, where $a$ denotes the axion field and $f_{a}$ is the axion
decay constant. Therefore, studies of the QCD interaction with axions at
finite chemical potential serve as studies of quark matter at finite $\theta$
and finite $\mu$ as well. Hence, within our work we will interchangeably use
$\theta$ in place of $a/f_{a}$.

Deep within the dense cores of compact stars, where matter is subjected to
extreme pressures and temperatures, exotic phases of strongly interacting
matter come into play, as the 2SC color superconductive phase [64, 65, 54, 66,
67]. This phase is characterized by a non-vanishing quark-quark condensate: as
such, it is not a color singlet, transforming as a color-antitriplet under
color rotations. In most previous calculations, in particular in those related
to the QCD phase diagram, only the scalar condensate was considered, as it is
the one that is favored by the one-instanton exchange. Within the present
work, we introduce both a scalar and a pseudo-scalar condensate, since both of
them are relevant when the coupling to the axions is considered. In agreement
with the common lore, we show that our model is consistent with favoring the
scalar condensate when the axion is not included in the model. On the other
hand, changing $\theta$ can result in a phase transition to a new ground state
where the condensate is a pseudo-scalar one, in qualitative agreement with
previous studies in normal quark matter [59, 49, 60, 61, 62]. Therefore, in
principle both condensates need to be introduced.

In the low-energy regime, where non-perturbative effects dominate, effective
models become indispensable tools for describing strongly interacting matter.
Chiral Perturbation Theory ($\chi$PT) stands out as a frequently employed
effective framework, significantly contributing to the understanding of the
vacuum structure of QCD and the properties of axions at low temperatures [68,
69, 70, 71, 72, 73]. $\chi$PT demonstrates notable advantages in low-energy
scenarios; for instance, its prediction of topological susceptibility at zero
temperature aligns well with lattice QCD findings [74, 75, 76]. However, its
applicability diminishes at high temperatures and/or large densities, as it
lacks information about QCD phase transitions. Hence, there’s a need for a
QCD-like model capable of accommodating axions and capturing the QCD phase
transition dynamics.One of the most popular is the Nambu-Jona-Lasinio (NJL)
model, which describes the dynamics of fermions and their interactions through
effective four-fermion interactions [53, 63]. The model includes the
instanton-induced interaction, which breaks the $U(1)_{A}$ symmetry, and it
can describe both the spontaneous breaking of chiral symmetry and how quarks
interact with axions.

The interaction of quark matter with axions has been already explored in [77,
59, 61, 49, 60, 62], in the framework of the NJL model, where the effect of
the chiral phase transition on the properties of the QCD axion was explored.
Furthermore, the axion potential was studied at finite quark chemical
potential, and the behavior of axion domain walls [78, 12, 79, 80] (see also
[81, 82] for a pedagogical introduction to walls) in bulk quark matter was
also investigated. It was found that the axion potential is very sensitive to
the chiral phase transition, particularly when quark matter is near
criticality, where the axion mass decreases and the self-coupling is enhanced.
In this study, we aim to extend our previous work [77] considering the
coupling of axions to diquarks, in order to take into account the effects of
the presence of axions in a color-superconducting medium. To our knowledge,
this is the first time that the coupling of axions to superconductive quarks
has been investigated in some detail. Our model is based on a four-fermion
effective interaction that contains a $U(1)_{A}-$preserving term, that can be
interpreted as an effective way to describe one-gluon exchange, as well as a
$U(1)_{A}-$breaking term, describing the instanton-mediated effective
interaction. The relative strength of the two terms is regulated by a
dimensionless, free parameter of the model, $\zeta$, while the overall
strength of the effective coupling, $G_{D}$, is chosen in order to reproduce a
given value of the superconductive gap when $\zeta$ is varied. Our main
results are related to the study of the full topological susceptibility,
$\chi=\partial^{2}\Omega/\partial\theta^{2}$, in the color-superconductive
phase: here, $\Omega$ denotes the full thermodynamic potential of
superconductive quark matter, while the derivative is understood at
$\theta=0$. We derive, within the model at hand, an exact analytical relation
between $\chi$ and the physical value of the superconductive gap. $\chi$ is
related to the axion mass, $m_{a}$, and decay constant, $f_{a}$, via
$\chi=m_{a}^{2}f_{a}^{2}$, therefore the knowledge of $\chi$ allows us to
extract an explicit, analytical formula for $m_{a}$ in the 2SC phase.
Moreover, we compute the full axion potential, and extract from it the low-
energy parameters of this potential: besides $m_{a}$, we compute the axion
self-coupling, $\lambda$. Within this work, we limit ourselves to a one-loop
approximation (usually called the mean field approximation), leaving the study
of the role of quantum fluctuations to future works.

The plan of the article is as follows. In section II we describe the model we
adopt in our work. In section III we present our results regarding the
superconductive gap at finite $\theta$ as well as the axion potential in the
two-flavor superconductive phase. Finally, in section IV we draw our
conclusions and discuss about possible future works. We use natural units
$\hbar=c=k_{B}=1$ throughout this article.

##  II The model

In this section, we present in some detail the model we use in our work. We
firstly present the Lagrangian density, describing how we couple the axion to
quarks in the superconductive phase of QCD. Next we turn to write the
thermodynamic potential, computed at one loop. We then show the shape of this
potential in the gaps space, discussing the location of its minima: this
discussion is useful to understand the results we show in the next section.
Finally, we present an approximated solution to the gap equation, valid for
small $\theta=a/f_{a}$, that allows us to qualitatively understand the
behavior of the superconductive gap versus $\theta$.

###  II.1 Lagrangian density

To begin with, we consider the Lagrangian density

$\displaystyle\mathcal{L}_{\mathrm{int}}$ $\displaystyle=$ $\displaystyle
g_{1}(q^{T}Ci\gamma_{5}\varepsilon\varepsilon
q)(\bar{q}i\gamma_{5}C\varepsilon\varepsilon\bar{q}^{T})+g_{2}(q^{T}C\varepsilon\varepsilon
q)(\bar{q}C\varepsilon\varepsilon\bar{q}^{T}).$ (1)

This has been used in many works on the superconductive phases of QCD, see for
example [58, 57]; in those references, it is assumed that $g_{1}=g_{2}$ so the
$U(1)_{A}$ symmetry is preserved at the Lagrangian level. In this work, we
assume $g_{1}\neq g_{2}$ from the very beginning, in order to break the
$U(1)_{A}$ symmetry at the level of the Lagrangian mimicking the instanton-
mediated interaction. In the Lagrangian density (1) we used
$C=-i\gamma_{2}\gamma_{0}$, satisfying $C^{2}=C^{\dagger}C=-1$. Moreover, we
adopted a condensed notation that suppresses the color and flavor indices
carried by the fields and the antisymmetric symbols in every bilinear. For
example,

$q^{T}Ci\gamma_{5}\varepsilon\varepsilon q=q^{T}_{\alpha
i}Ci\gamma_{5}\varepsilon_{\alpha\beta 3}\varepsilon_{ij3}q_{\beta j};$ (2)

here $\alpha$, $\beta$ denote color indices, while $i$, $j$ stand for flavor
indices.

Starting from (1), we isolate a term that is invariant under $U(1)_{A}$ and a
term that explicitly breaks this symmetry. This can be easily achieved by
adding and subtracting the terms
$g_{2}(q^{T}Ci\gamma_{5}\varepsilon\varepsilon
q)(\bar{q}i\gamma_{5}C\varepsilon\varepsilon\bar{q}^{T})$ and
$g_{1}(q^{T}C\varepsilon\varepsilon
q)(\bar{q}C\varepsilon\varepsilon\bar{q}^{T})$; we then get

$\displaystyle\mathcal{L}_{\mathrm{int}}$ $\displaystyle=$
$\displaystyle\frac{(g_{1}+g_{2})}{2}\left[(q^{T}Ci\gamma_{5}\varepsilon\varepsilon
q)(\bar{q}i\gamma_{5}C\varepsilon\varepsilon\bar{q}^{T})+(q^{T}C\varepsilon\varepsilon
q)(\bar{q}C\varepsilon\varepsilon\bar{q}^{T})\right]$ (3) $\displaystyle+$
$\displaystyle\frac{(g_{1}-g_{2})}{2}\left[(q^{T}Ci\gamma_{5}\varepsilon\varepsilon
q)(\bar{q}i\gamma_{5}C\varepsilon\varepsilon\bar{q}^{T})-(q^{T}C\varepsilon\varepsilon
q)(\bar{q}C\varepsilon\varepsilon\bar{q}^{T})\right].$

The first line in the above Lagrangian density is now $U(1)_{A}-$preserving,
so we confined the breaking of $U(1)_{A}$ to the second line of (3). Writing
$q=(\mathcal{P}_{L}+\mathcal{P}_{R})q$, where

$\mathcal{P}_{R}=\frac{1+\gamma_{5}}{2},~{}~{}~{}\mathcal{P}_{L}=\frac{1-\gamma_{5}}{2},$
(4)

we can easily rewrite the second line of (3) as

$(g_{1}-g_{2})\left[(q^{T}Ci\gamma_{5}\mathcal{P}_{L}\varepsilon\varepsilon
q)(\bar{q}i\gamma_{5}\mathcal{P}_{L}C\varepsilon\varepsilon\bar{q}^{T})+(q^{T}Ci\gamma_{5}\mathcal{P}_{R}\varepsilon\varepsilon
q)(\bar{q}i\gamma_{5}\mathcal{P}_{R}C\varepsilon\varepsilon\bar{q}^{T})\right],$
(5)

where we used $\gamma_{5}\mathcal{P}_{L}q=-\mathcal{P}_{L}q=-q_{L}$ and
$\gamma_{5}\mathcal{P}_{R}q=\mathcal{P}_{R}q=q_{R}$. Finally, defining
$G_{D}=(g_{1}+g_{2})/2$ and $\zeta G_{D}=(g_{1}-g_{2})$, we can rewrite the
interaction (1) as

$\displaystyle\mathcal{L}_{\mathrm{int}}$ $\displaystyle=$ $\displaystyle
G_{D}\left[(q^{T}Ci\gamma_{5}\varepsilon\varepsilon
q)(\bar{q}i\gamma_{5}C\varepsilon\varepsilon\bar{q}^{T})+(q^{T}C\varepsilon\varepsilon
q)(\bar{q}C\varepsilon\varepsilon\bar{q}^{T})\right]$ (6) $\displaystyle+\zeta
G_{D}\left[(q^{T}Ci\gamma_{5}\mathcal{P}_{L}\varepsilon\varepsilon
q)(\bar{q}i\gamma_{5}\mathcal{P}_{L}C\varepsilon\varepsilon\bar{q}^{T})+(q^{T}Ci\gamma_{5}\mathcal{P}_{R}\varepsilon\varepsilon
q)(\bar{q}i\gamma_{5}\mathcal{P}_{R}C\varepsilon\varepsilon\bar{q}^{T})\right],$

The second line in (6) breaks $U(1)_{A}$: we consider it as an effective way
to model the quark-quark interaction arising from the one-instanton exchange.

In order to couple the axions to the quarks, we note that the former can only
couple to the $U(1)_{A}$-breaking term in Eq. (6) (in fact, the QCD axion
couples to instanton-like gluon configurations). In agreement with what has
been done for the coupling to quarks in the vacuum [77], we write this
coupling as

$\displaystyle\mathcal{L}_{\mathrm{int}}$ $\displaystyle=$ $\displaystyle
G_{D}\left[(q^{T}Ci\gamma_{5}\varepsilon\varepsilon
q)(\bar{q}i\gamma_{5}C\varepsilon\varepsilon\bar{q}^{T})+(q^{T}C\varepsilon\varepsilon
q)(\bar{q}C\varepsilon\varepsilon\bar{q}^{T})\right]$ (7) $\displaystyle+\zeta
G_{D}\left[e^{i\frac{a}{f_{a}}}(q^{T}Ci\gamma_{5}\mathcal{P}_{L}\varepsilon\varepsilon
q)(\bar{q}i\gamma_{5}\mathcal{P}_{L}C\varepsilon\varepsilon\bar{q}^{T})+e^{-i\frac{a}{f_{a}}}(q^{T}Ci\gamma_{5}\mathcal{P}_{R}\varepsilon\varepsilon
q)(\bar{q}i\gamma_{5}\mathcal{P}_{R}C\varepsilon\varepsilon\bar{q}^{T})\right].$

The Lagrangian density (7) specifies the interaction we adopt in our model. As
emphasized above, the first line in the right hand side of (7) preserves the
$U(1)_{A}$ symmetry, representing an effective way of writing the attractive
channel of the quark-quark interaction arising from one-gluon exchange. The
axial symmetry instead is broken by the terms in the second line [83, 84]. We
note that in [58], authors included both terms in the first line of (7);
however, they usually neglect the second addendum in the first line of (7),
since in the mean-field, it gives rise to a pseudo-scalar condensate that is
usually neglected: as a matter of fact the instanton-induced interaction,
namely the second line of (7), favors condensation in the scalar channel when
$a=0$ (we explicitly verified this statement within our model, see below).
However, when $a\neq 0$ condensation can happen in the pseudo-scalar channel
as well, therefore we need to consider the whole interaction (7).

In this work, we treat the interaction (7) within the mean field
approximation. In order to implement this, we introduce the condensates

$\displaystyle\langle q^{T}Ci\gamma_{5}\mathcal{P}_{L}\varepsilon\varepsilon
q\rangle=-h_{L},$ (8)
$\displaystyle\langle\bar{q}i\gamma_{5}\mathcal{P}_{L}C\varepsilon\varepsilon\bar{q}^{T}\rangle=h_{R}^{*},$
(9) $\displaystyle\langle
q^{T}Ci\gamma_{5}\mathcal{P}_{R}\varepsilon\varepsilon q\rangle=h_{R},$ (10)
$\displaystyle\langle\bar{q}i\gamma_{5}\mathcal{P}_{R}C\varepsilon\varepsilon\bar{q}^{T}\rangle=-h_{L}^{*},$
(11)

as well as their combinations

$\displaystyle\langle q^{T}Ci\gamma_{5}\varepsilon\varepsilon
q\rangle=h_{R}-h_{L},$ (12)
$\displaystyle\langle\bar{q}i\gamma_{5}C\varepsilon\varepsilon\bar{q}^{T}\rangle=h_{R}^{*}-h_{L}^{*}.$
(13)

Using these, as well as the assumptions $h_{L}=h_{L}^{*}$, $h_{R}=h_{R}^{*}$,
we get, within the mean-field approximation,

$\displaystyle\mathcal{L}_{\mathrm{int}}$ $\displaystyle=$ $\displaystyle
G_{D}\varepsilon\varepsilon(h_{R}-h_{L})\left[(\bar{q}i\gamma_{5}C\varepsilon\varepsilon\bar{q}^{T})+(q^{T}Ci\gamma_{5}\varepsilon\varepsilon
q)\right]-G_{D}(h_{R}-h_{L})^{2}$ (14) $\displaystyle-
G_{D}\varepsilon\varepsilon\left[(h_{R}+h_{L})(\bar{q}iC\varepsilon\varepsilon\bar{q}^{T})-(h_{R}+h_{L})(q^{T}C\varepsilon\varepsilon
q)\right]-G_{D}(h_{R}+h_{L})^{2}$ $\displaystyle+\zeta
G_{D}\varepsilon\varepsilon
e^{i\frac{a}{f_{a}}}\left[-h_{L}(\bar{q}i\gamma_{5}\mathcal{P}_{L}C\varepsilon\varepsilon\bar{q}^{T})+h_{R}(q^{T}Ci\gamma_{5}\mathcal{P}_{L}\varepsilon\varepsilon
q)\right]+\zeta G_{D}e^{i\frac{a}{f_{a}}}h_{L}h_{R}$ $\displaystyle+\zeta
G_{D}\varepsilon\varepsilon
e^{-i\frac{a}{f_{a}}}\left[h_{R}(\bar{q}i\gamma_{5}\mathcal{P}_{R}C\varepsilon\varepsilon\bar{q}^{T})-h_{L}(q^{T}Ci\gamma_{5}\mathcal{P}_{R}\varepsilon\varepsilon
q)\right]+\zeta G_{D}e^{-i\frac{a}{f_{a}}}h_{L}h_{R}.$

Note that $\mathcal{L}_{\mathrm{int}}=\mathcal{L}_{\mathrm{int}}^{\dagger}$.

In order to simplify the notation we introduce the Nambu-Gorkov bi-spinors

$\Psi=\left(\begin{array}[]{c}q\\\
C\bar{q}^{T}\end{array}\right),\qquad\bar{\Psi}=\left(\bar{q},q^{T}C\right).$
(15)

In terms of these, we then can rewrite the interaction lagrangian density as

$\mathcal{L}_{\mathrm{int}}=\mathcal{V}+\bar{\Psi}\bm{\Delta}\Psi,$ (16)

where

$\mathcal{V}=-2G_{D}(h_{R}^{2}+h_{L}^{2})+2\zeta
G_{D}h_{L}h_{R}\cos\left(\frac{a}{f_{a}}\right),$ (17)

and

$\bm{\Delta}=\left(\begin{array}[]{cc}0&\Phi^{-}\\\ &\\\
\Phi^{+}&0\end{array}\right),$ (18)

with

$\displaystyle\Phi^{-}$ $\displaystyle=$ $\displaystyle
G_{D}\left[2h_{R}\mathcal{P}_{L}-2h_{L}\mathcal{P}_{R}-e^{ia/f_{a}}h_{L}\zeta\mathcal{P}_{L}+e^{-ia/f_{a}}h_{R}\zeta\mathcal{P}_{R}\right]i\gamma_{5}\varepsilon_{ij}\varepsilon_{\alpha\beta
3}$ (19) $\displaystyle\Phi^{+}$ $\displaystyle=$ $\displaystyle
G_{D}\left[2h_{R}\mathcal{P}_{R}-2h_{L}\mathcal{P}_{L}+e^{ia/f_{a}}h_{R}\zeta\mathcal{P}_{L}-e^{-ia/f_{a}}h_{L}\zeta\mathcal{P}_{R}\right]i\gamma_{5}\varepsilon_{ij}\varepsilon_{\alpha\beta
3}.$ (20)

In the above equations, the first two terms arise from the one-gluon exchange
interaction, while those proportional to $\zeta$ arise from the one-instanton
exchange. Moreover, $\bm{1}_{D}$ denotes the identity matrix in the Dirac
space. The matrix $\bm{\Delta}$ has a similar structure to that of [57]:
indeed, it satisfies $\Phi^{+}=\gamma_{0}(\Phi^{-})^{\dagger}\gamma_{0}$.

To the interaction (16) we need to add the kinetic term of quarks at finite
chemical potential $\mu$: this contribution is well-known [58, 57] and leads
to the full lagrangian, that in momentum space reads

$\mathcal{L}=\mathcal{V}+\bar{\Psi}S^{-1}\Psi.$ (21)

Here, the inverse quark propagator is given by

$S^{-1}(p)=\left(\begin{array}[]{cc}(\not{p}+\mu\gamma_{0})\bm{1}_{C}\bm{1}_{F}&\Phi^{-}\\\
&\\\ \Phi^{+}&(\not{p}-\mu\gamma_{0})\bm{1}_{C}\bm{1}_{F}\end{array}\right),$
(22)

where $\bm{1}_{C}$ and $\bm{1}_{F}$ correspond to the identities in color and
flavor spaces respectively. We note that these matrices have dimension
$4\times 3\times 2\times 2=48$, due to Dirac, color, flavor, and Gorkov
indices, respectively.

###  II.2 Thermodynamic potential

The thermodynamic potential is obtained via the standard integration over the
fermion fields in the partition function [58, 57], which leads at

$\Omega=-\mathcal{V}+\Omega_{1-\mathrm{loop}},$ (23)

where $\mathcal{V}$ denotes the mean field contribution (17), and
$\Omega_{1-\mathrm{loop}}$ corresponds to the 1-loop contribution of the
quarks, namely

$\Omega_{1-\mathrm{loop}}=-T\sum_{n}\int\frac{d^{3}\,p}{(2\pi)^{3}}\frac{1}{2}\mathrm{Tr}\log\left(\beta
S^{-1}(i\omega_{n},\vec{p})\right).$ (24)

In getting Eq. (24) we performed the functional integral over fermions at
fixed $a$: in fact, in this work the axion field is treated as a classical
background. Consequently, the thermodynamic potential is a function of $a$ and
physical quantities have to be computed at a given value of $a$. Moreover, we
used the imaginary time formalism of finite temperature field theory, and
$\omega_{n}=(2n+1)\pi T$ are the relevant fermionic Matsubara frequencies. The
overall $1/2$ in (24) takes into account the artificial doubling of degrees of
freedom introduced by shifting to the Nambu-Gorkov bi-spinors. Finally, the
trace is understood in Nambu-Gorkov, Dirac, color and flavor spaces.

In order to evaluate the sum over $\omega_{n}$ in Eq. (24), we follow the
well-known strategy of computing the eigenvalues of the matrix [58, 57]

$\mathcal{T}=\left(\begin{array}[]{cc}(-\gamma_{0}\vec{p}\cdot\vec{\gamma}+\mu)\bm{1}_{C}\bm{1}_{F}&\gamma_{0}\Phi^{-}\\\
&\\\
\gamma_{0}\Phi^{+}&(-\gamma_{0}\vec{p}\cdot\vec{\gamma}-\mu)\bm{1}_{C}\bm{1}_{F}\end{array}\right).$
(25)

These are given by

$\displaystyle\varepsilon_{1,\pm}$ $\displaystyle=$ $\displaystyle\pm|p-\mu|,$
(26) $\displaystyle\varepsilon_{2,\pm}$ $\displaystyle=$
$\displaystyle\pm|p+\mu|,$ (27) $\displaystyle\varepsilon_{3,\pm}$
$\displaystyle=$ $\displaystyle\pm\sqrt{(p-\mu)^{2}+\Delta_{3}^{2}},$ (28)
$\displaystyle\varepsilon_{4,\pm}$ $\displaystyle=$
$\displaystyle\pm\sqrt{(p+\mu)^{2}+\Delta_{3}^{2}},$ (29)
$\displaystyle\varepsilon_{5,\pm}$ $\displaystyle=$
$\displaystyle\pm\sqrt{(p-\mu)^{2}+\Delta_{5}^{2}},$ (30)
$\displaystyle\varepsilon_{6,\pm}$ $\displaystyle=$
$\displaystyle\pm\sqrt{(p+\mu)^{2}+\Delta_{5}^{2}},$ (31)

where $p\equiv|\bm{p}|$ denotes the magnitude of the $3-$momentum $\bm{p}$;
each of the eigenvalues above has a multiplicity equal to $4$. The eigenvalues
(26) and (27) correspond to those of the blue quarks, which do not participate
to the pairing and remain ungapped. Moreover, we put

$\displaystyle\Delta_{3}^{2}$ $\displaystyle=$
$\displaystyle\zeta^{2}G_{D}^{2}h_{L}^{2}+4G_{D}^{2}h_{R}^{2}-4\zeta
G_{D}^{2}h_{L}h_{R}\cos(a/fa),$ (32) $\displaystyle\Delta_{5}^{2}$
$\displaystyle=$
$\displaystyle\zeta^{2}G_{D}^{2}h_{R}^{2}+4G_{D}^{2}h_{L}^{2}-4\zeta
G_{D}^{2}h_{L}h_{R}\cos(a/fa),$ (33)

that represent the (squared) gaps in the quark spectrum.

Putting $\varepsilon_{k}(\vec{p})=\varepsilon_{k,+}(\vec{p})$, $k=1,\cdots,6$,
we use $\ln\det\beta S^{-1}=\mbox{Tr}\ln\beta S^{-1}$ and write

$\ln\det\beta
S^{-1}_{L}(i\omega_{n},\vec{p})=4\sum_{k=1}^{6}\ln\left(\frac{\omega_{n}^{2}+\varepsilon_{k}(\vec{p})^{2}}{T^{2}}\right),$
(34)

where we made explicit the degeneracy, equal to $4$, of the eigenvalues and
matched the ones with opposite signs. Now the Matsubara sum can be easily
evaluated, since

$T\sum_{n}\ln\left(\frac{\omega_{n}^{2}+\varepsilon_{k}(\vec{p})^{2}}{T^{2}}\right)=|\varepsilon_{k}(\vec{p})|+2T\ln\left(1+e^{-|\varepsilon_{k}(\vec{p})|/T}\right),$
(35)

so the 1-loop contribution finally results in

$\Omega_{1-\mathrm{loop}}^{L}=-2\sum_{k=1}^{6}\int\frac{d^{3}p}{(2\pi)^{3}}\left[\varepsilon_{k}(\vec{p})+2T\ln\left(1+e^{-\varepsilon_{k}(\vec{p})/T}\right)\right],$
(36)

where we took into account that according to our definitions the
$\varepsilon_{k}$ are always positive. Therefore, the thermodynamic potential
is given by

$\Omega=2G_{D}(h_{R}^{2}+h_{L}^{2})-2\zeta
G_{D}h_{L}h_{R}\cos\left(\frac{a}{f_{a}}\right)-2\sum_{k=1}^{6}\int\frac{d^{3}p}{(2\pi)^{3}}\left[\varepsilon_{k}(\vec{p})+2T\ln\left(1+e^{-\varepsilon_{k}(\vec{p})/T}\right)\right].$
(37)

The first addendum in the right-hand side of (36) corresponds to the $T=0$
contribution to the thermodynamic potential, which is divergent in the
ultraviolet. In order to regulate this divergence, we introduce a sharp 3D
momentum cutoff, $\Lambda$, and rewrite that term as

$\int\frac{d^{3}p}{(2\pi)^{3}}\varepsilon_{k}(p)=\frac{4\pi}{8\pi^{3}}\int_{0}^{\Lambda}p^{2}dp~{}\varepsilon_{k}(p).$
(38)

The cutoff roughly represents the momentum scale above which the contact
interaction used in the present work should replaced by a non-local term
directly borrowed from QCD. We treat $\Lambda$ as a free parameter of the
model, and from previous studies based on the NJL model, we assume
$\Lambda=O(1~{}\mathrm{GeV})$.

We notice that although in (36) we sum over the six positive eigenvalues of
$\mathcal{T}$, the ones corresponding to the blue quarks, which do not
participate in the pairing, do not contribute to the value of the condensates.

It is useful to stress that both $\mathcal{V}$ and the quark spectrum are
invariant under the set of transformations $h_{L}\leftrightarrow h_{R}$ and
$h_{L}\leftrightarrow-h_{R}$; moreover, if $\zeta=0$ then $\Omega$ depends on
$h_{L}^{2}$ and $h_{R}^{2}$ only: as a consequence, in this limit we expect
$\Omega$ to develop a set of four degenerate minima along the lines $h_{L}=\pm
h_{R}$. This degeneracy is removed by $\zeta\neq 0$. This picture is confirmed
by the direct evaluation of $\Omega$, see Fig. 1 in the next section.

###  II.3 Thermodynamic potential at finite $\theta=a/f_{a}$: the lines
$h_{L}=\pm h_{R}$

Figure 1: Thermodynamic potential at $T=0$, computed for $\theta\equiv
a/f_{a}=0$, and for $\zeta=0$ (left panel) and for $\zeta=0.2$ (right panel)
In the plots, $x$ and $y$-axes denote $\Delta_{L}$ and $\Delta_{R}$
respectively (measured in MeV).

In this subsection, we analyze the shape of $\Omega$ in the $(h_{L},h_{R})$
plane: this is preparatory to the results on the gap that we present in the
next section. In Fig. 1 we plot the thermodynamic potential at $T=0$, computed
for $\theta\equiv a/f_{a}=0$, and for $\zeta=0$ (left panel) and for
$\zeta=0.2$ (right panel). In the plots, $x$ and $y$-axes denote $\Delta_{L}$
and $\Delta_{R}$ respectively (measured in MeV), which are defined as

$\Delta_{L}=2G_{D}h_{L},~{}~{}~{}\Delta_{R}=2G_{D}h_{R}.$ (39)

It is more convenient to use $\Delta_{L,R}$ instead of $h_{L,R}$ since the
former correspond to the gaps in the quark spectrum when $\zeta=0$. From
these, the scalar and the pseudoscalar condensate can be formed as

$\Delta_{\mathrm{S}}=\Delta_{R}-\Delta_{L},~{}~{}~{}\Delta_{\mathrm{PS}}=\Delta_{R}+\Delta_{L}.$
(40)

We notice in Fig. 1 that for $\zeta=0$, namely when we consider the one-gluon-
exchange interaction only, there are two loci of degenerate minima, namely
$\Delta_{L}=\Delta_{R}$ that corresponds to the pseudoscalar condensate and
$\Delta_{L}=-\Delta_{R}$ corresponding to the scalar condensate. On the other
hand, when the instanton-induced interaction is switched on $(\zeta\neq 0)$
this degeneracy is removed, and the condensation in the scalar channel
$\Delta_{L}=-\Delta_{R}$ is favored.

Figure 2: Thermodynamic potential at $T=0$, computed for $a_{1}=a_{2}=1$,
$\zeta=0.2$, and several values of $\theta\equiv a/f_{a}$. In each plot, $x$
and $y$-axes denote $\Delta_{L}$ and $\Delta_{R}$ respectively (measured in
MeV). Top-left plot corresponds to $\theta=0$, top-right to
$\theta=\pi/2-\varepsilon$, center-left to $\theta=\pi/2$, center-right to
$\theta=\pi/2+\varepsilon$, bottom-left to $\theta=\pi$, and finally bottom-
right to $\theta=2\pi$.

In Fig. 2 we plot the thermodynamic potential at $T=0$, computed for
$\zeta=0.2$, and several values of $\theta\equiv a/f_{a}$. In each plot, $x$
and $y$-axes denote $\Delta_{L}$ and $\Delta_{R}$ respectively (measured in
MeV). In the figure, the top-left plot corresponds to $\theta=0$, top-right to
$\theta=\pi/2-\varepsilon$, center-left to $\theta=\pi/2$, center-right to
$\theta=\pi/2+\varepsilon$, bottom-left to $\theta=3\pi/2$, and finally
bottom-right to $\theta=2\pi$. For all the cases shown in the figure, the
minima of $\Omega$ sit on the lines $\Delta_{L}=\pm\Delta_{R}$: this is the
obvious consequence of the fact that $\Omega$ is invariant under the
transformations $h_{L}\leftrightarrow\pm h_{R}$. For $a/f_{a}$ in the range
$[0,\pi/2)$, the global minima correspond to $\Delta_{L}=-\Delta_{R}$; hence,
in this range of $a/f_{a}$ only the scalar condensate exists. The minima are
dried out as $a/f_{a}$ is increased: indeed, for $a/f_{a}$ in the range
$(\pi/2,3\pi/2)$ the global minima are located along the line
$\Delta_{L}=\Delta_{R}$, implying a phase transition from the scalar to the
pseudo-scalar condensate. Finally, increasing $a/f_{a}$ up to $2\pi$ the
location of the minima changes again, resulting in another phase transition
from the pseudo-scalar to the scalar condensate. Along the lines
$\Delta_{L}=\pm\Delta_{R}$ we have $\Delta_{3}=\Delta_{5}$ for any $a$ and
$\zeta$: hence, there is only one gap in the quark spectrum, see (32) and
(33).

We notice that the results on the minima of $\Omega$ discussed above stand for
$\zeta>0$: we checked that for $\zeta<0$ the role of the scalar and the
pseudoscalar condensates invert; besides this, there is no major difference
between the system with positive and negative $\zeta$. Consequently, from now
on we limit ourselves to show results for $\zeta>0$ only.

Before going ahead, we determine the allowed values of $\zeta$ that lead to
the nontrivial solution of the gap equation. To this end, it is enough to
limit ourselves to the gap equation at $T=0$ and small $\theta$. Imposing the
stationarity condition of $\Omega$ with respect to $h_{L}$, namely
$\partial\Omega/\partial h_{L}=0$ with $\Omega$ given by Eq. (37), and then
specializing this condition to the line $h_{L}=-h_{R}$, which is the one along
which $\Omega$ develops minima for small $\theta$ as discussed above, gives

$2+\zeta\cos(a/f_{a})=\frac{G_{D}}{2\pi^{2}}\int_{0}^{\Lambda}p^{2}dp~{}\left(\frac{A(a,\zeta)}{\sqrt{(p-\mu)^{2}+G_{D}^{2}h_{L}^{2}A(a,\zeta)}}+\mu\rightarrow-\mu\right),$
(41)

where we removed the trivial solution $h_{L}=0$, and we defined

$A(a,\zeta)=4+4\zeta\cos(a/f_{a})+\zeta^{2}.$ (42)

The above equation shows that the condition $\zeta<2$ must be satisfied in
order to have a non-trivial solution of the gap equation. In fact, the right-
hand side of (41) is always positive, and so the left-hand side must be
positive as well: this can be obtained if and only if $\zeta<2$. For
$\zeta\geq 2$ there is a range of $a$ in which the gap equation has only the
trivial solution $h_{L}=0$: hence, in what follows we limit ourselves to
$\zeta<2$.

##  III Results

###  III.1 Gaps versus $a/f_{a}$

Figure 3: Gaps $\Delta_{L}$ and $\Delta_{R}$ at $T=0$ and $\mu=400$ MeV,
versus $a/f_{a}$, and for several values of $\zeta$. $G_{D}$ is fixed in order
to have $\Delta_{L}=50$ MeV for $a=0$ at $\mu=400$ MeV.

In this subsection, we present our results for the gap parameters. For
$\Delta_{L}=\pm\Delta_{R}$, which correspond to the directions along which the
minima develop at finite $a$, we get from Eqs. (32) and (33) that
$\Delta_{3}=\Delta_{5}$. Therefore we limit ourselves to show results for
$\Delta_{3}$ only. We fix $\mu=400$ MeV and $\Lambda=1$ GeV as representative
values of the quark chemical potential and the UV cutoff. For each $\zeta$, we
fix $G_{D}$ in order to have a desired value of $\Delta_{L}$ at $a=0$.

In Fig. 3 we plot $\Delta_{L}$ and $\Delta_{R}$ versus $a/f_{a}$ at $T=0$ and
$\mu=400$ MeV; we consider several values of $\zeta$, while $G_{D}$ is fixed
for each $\zeta$ so that $\Delta_{L}=50$ MeV at $a=0$. We firstly focus on the
results for $\zeta=0.2$, $0.5$ and $1$. For these values of $\zeta$, we notice
that for $a/f_{a}$ in the range $(0,\pi/2)$, $\Delta_{R}=-\Delta_{L}$, hence
only the scalar condensate forms. On the other hand, for $a/f_{a}$ in the
range $(\pi/2,3\pi/2)$, $\Delta_{R}=\Delta_{L}$: in this case, only the
pseudo-scalar condensate forms. Finally, for $a/f_{a}$ in the range
$(3\pi/2,2\pi)$ we find again $\Delta_{R}=-\Delta_{L}$, hence there is a phase
transition from the pseudo-scalar to the scalar condensate. We also notice
that for $\zeta=1.75$ the magnitude of $\Delta_{L}$ increases with $a$ for
$a/f_{a}$ in the range $(0,\pi/2)$.

The behavior of $\Delta_{L}$ near $a=0$ can be semi-quantitatively understood
by the approximate solution of the gap equation (41), namely

$\Delta_{L}=\frac{4\delta}{\sqrt{A(a,\zeta)}}\exp\left(-\frac{\pi^{2}}{\mu^{2}G_{D}\Theta(a,\zeta)}\right),$
(43)

where $A$ is given by Eq. (42) and

$\Theta(a,\zeta)=\frac{A(a,\zeta)}{2+\zeta\cos(a/f_{a})}.$ (44)

The approximated solution (43) can be obtained within the High Density
Effective Theory (HDET) of QCD, see [85] and references therein, as well as
appendix A; it has the standard form of the Bardeen-Cooper-Schrieffer (BCS)
gap in the theory of superconductivity [86, 87]: in fact, $\mu^{2}$ is
proportional to the density of states of the pairing quarks at the Fermi
surface, and $\delta$ in (43) plays the role of the Debye frequency,
$\omega_{D}$, of the BCS theory, which cuts high momentum modes out of the
pairing. The solution (43) is strictly valid only in the weak coupling, hence
it cannot quantitatively reproduce the results in Fig. 3: however, it is still
helpful to graps the behavior of $\Delta_{L}$ near $a=0$. In fact, from (43)
we notice that for a fixed value of $\zeta$, the coupling of the
superconductive quarks to the axion field effectively affects the pairing in
two ways. Firstly, it changes the width of the shell around the Fermi surface
that contributes to the pairing, namely

$\frac{2\delta}{\sqrt{A(0,\zeta)}}\rightarrow\frac{2\delta}{\sqrt{A(a,\zeta)}}.$
(45)

Secondly, it effectively changes the chemical potential of the quarks, that is

$\mu^{2}\Theta(0,\zeta)\rightarrow\mu^{2}\Theta(a,\zeta).$ (46)

Both $A(a,\zeta)$ and $\Theta(a,\zeta)$ are decreasing functions of $a$ at
fixed $\zeta$. Consequently, the response of $\Delta_{L}$ to the coupling with
the axion has two competing effects: on the one hand, it leads to the opening
of the shell around the Fermi surface, thus increasing the portion of phase
space involved in the pairing; on the other hand, it effectively lowers the
chemical potential of the paired quarks, implying the decrease of the volume
of phase space available for pairing. Which of the two effects wins depends,
at a given $G_{D}$, on the value of $\zeta$. This can be seen by inspecting
the curvature of $\Delta_{L}$ for small $\theta=a/f_{a}$. To this end, it is
enough to expand Eq. (43) near $\theta=0$, getting

$\Delta_{L}=\Delta_{L,0}\left(1+\frac{\kappa}{2}\theta^{2}\right),$ (47)

where $\Delta_{L,0}$ denotes the gap (43) at $\theta=0$, and the curvature is

$\kappa=\zeta\frac{\pi^{2}(\zeta-2)+2G_{D}\mu^{2}(2+\zeta)}{G_{D}\mu^{2}(2+\zeta)^{4}}.$
(48)

The curvature is trivially zero for $\zeta=0$, as well as for

$\zeta=\zeta_{\mathrm{crit}}=\frac{2(\pi^{2}-G_{D}\mu^{2})}{\pi^{2}+G_{D}\mu^{2}};$
(49)

it is negative for $\zeta$ in the range $(0,\zeta_{\mathrm{crit}})$, positive
otherwise. Hence, the response of $\Delta_{L}$ to $\theta$ around $\theta=0$
depends on the specific value of $\zeta$, and the turning point of $\zeta$,
namely $\zeta_{\mathrm{crit}}$, depends on the value of $G_{D}$.

Figure 4: Gap in the quark spectrum, $\Delta_{3}$, at $T=0$ and $\mu=400$ MeV,
versus $a/f_{a}$, and for several values of $\zeta$. $G_{D}$ is fixed in order
to have $\Delta_{L}=50$ MeV for $a=0$ at $\mu=400$ MeV.

The gap in the quark spectrum, along the global minima lines, is given by
$\Delta_{3}$ in (32). We show $\Delta_{3}$ versus $a/f_{a}$ in Fig. 4. We
notice that $\Delta_{3}$ is a periodic function of $a/f_{a}$, with a period
equal to $\pi$, as expected from the results on $\Delta_{L}$ and $\Delta_{R}$
shown Fig. 3. We find that although the response of $\Delta_{L}$ and
$\Delta_{R}$ on $a/f_{a}$ depends on $\zeta$, hence on the weight of the
t’Hooft interaction term, the qualitative behavior of $\Delta_{3}$ does not
depend on $\zeta$.

###  III.2 Axion potential

Figure 5: Axion potential, $V(\theta)=\Omega(\theta)-\Omega(0)$ with
$\theta=a/f_{a}$, versus $a/f_{a}$, computed at $T=0$ and $\mu=400$ MeV and
for several values of $\zeta$. The potential is measured in units of
$\Lambda^{4}$ with $\Lambda=1$ GeV. $G_{D}$ is fixed to have $\Delta_{L}=50$
MeV at $a=0$.

One of the main results of our work is the computation of the axion potential,
$V(a)$, in a superconductive phase of QCD. In Fig. 5 we plot the
$V(a)\equiv\Omega(a)-\Omega(a=0)$ versus $a/f_{a}$, at $T=0$ and $\mu=400$
MeV. For each value of $a/f_{a}$, the potential has been computed at the
global minimum in the $(\Delta_{L},\Delta_{R})$ space. The behavior of $V$ is
in agreement with that of the quark spectrum shown in Fig. 4, particularly for
what concerns the periodicity. We notice that $a=0$ is a global minimum of
$V(a)$ for the whole range of parameters studied here. Other degenerate minima
are located at $a/f_{a}=n\pi/2$ with $n=\pm 1,\pm 2,\dots$.

It is useful to stress the difference with respect to the axion potential
computed in presence of the chiral and the $\eta-$condensate [77]: indeed, in
that case, the periodicity of $V(a)$ is equal to $2\pi$, while in the case of
superconductive matter, we find a periodicity equal to $\pi$. This can be
understood easily. In fact, in the present case at $a=0$ only the scalar
condensate forms, then increasing the value of $a$ the global minima dry out
and at $a/f_{a}=\pi/2$ the thermodynamic potential develops four degenerate
minima, see Fig. 2; further increasing $a/f_{a}$ in the range $(\pi/2,3\pi/2)$
results in the rotation of the global minima of $\Omega$ and in the consequent
formation of the pseudo-scalar condensate. The ground states at $a=0$ and
$a/f_{a}=\pi$ are degenerate, hence the periodicity of the potential, but the
two adjacent minima correspond to two different condensation channels: in
particular, the ground state at $a=0$ is characterized by scalar condensation,
while the minimum at $\theta=\pi/2$ corresponds to condensation in the pseudo-
scalar channel.

###  III.3 Topological susceptibility

Figure 6: Topological susceptibility, $\chi^{1/4}$, versus $\zeta$, at $T=0$
and for several values of $\Delta_{L}$. We used $\mu=400$ MeV, $\Lambda=1000$
MeV. The value of $\Delta_{L}$ in the legend has been used to fix the value of
$G_{D}$, so that for each line $\Delta_{L}$ is kept fixed while $\zeta$ is
changed. The analytical curve corresponds to Eq. (52) computed for
$\Delta_{L}=50$ MeV and $G_{D}\Lambda^{2}=4.69$.

In this section, we analyze the full topological susceptibility, $\chi$, which
measures the fluctuations of the topological charge and is defined as

$\chi=\left.\frac{d^{2}\Omega}{d\theta^{2}}\right|_{\theta=0},~{}~{}~{}\theta=\frac{a}{f_{a}}.$
(50)

In the numerical calculation, we treat the above derivative as a total
derivative, that in principle takes into account also of the dependence of the
condensates on $\theta$ at $\theta=0$.

In Fig. 6 we plot the fourth root of versus $\chi$, at $T=0$ and for several
values of $\Delta_{L}$. We used $\mu=400$ MeV, $\Lambda=1000$ MeV. The value
of $\Delta_{L}$ in the legend has been used to fix the value of $G_{D}$, so
that for each line $\Delta_{L}$ is kept fixed while $\zeta$ is changed. $\chi$
is an interesting quantity by itself, since it encodes information about the
fluctuations of the topological charge in the dense and superconductive QCD
medium. Moreover, it is directly related to the squared axion mass, see the
next subsection. Previous studies performed within Lattice QCD in the isospin-
symmetric case, as well as within chiral perturbation theory and NJL models,
agree on the value $\chi^{1/4}\approx 78$ MeV at $T=0$ and $\mu=0$ [74, 8, 59,
88], as well as with the Di Vecchia-Crewther-Leutwyler-Smilga-Veneziano
formula [89, 90, 91, 92], that for two flavors reads

$\chi=|\langle\bar{q}q\rangle|\frac{m_{u}m_{d}}{m_{u}+m_{d}},$ (51)

where $\langle\bar{q}q\rangle$ is the chiral condensate in the vacuum. $\chi$
then decreases at high $T$ where the smooth crossover to the quark-gluon
plasma phase takes place [74, 8, 59, 88].

The results in Fig. 6 show that for a given value of $\zeta$, the topological
susceptibility increases with the strength of the coupling, as expected.
Moreover, if the coupling is tuned in order to give a value of $\Delta_{L}$,
changing $\zeta$ within the range $(0.5,1.5)$ does not substantially affect
$\chi^{1/4}$. We also note that keeping $\zeta$ in the aforementioned range
keeps $\chi^{1/4}$ in the superconductive phase in the same ballpark of the
value it takes in the vacuum, unless we take a very large superconductive gap
as in the case $\Delta_{L}=100$ MeV shown in Fig. 6.

Within our model we were able to obtain an analytical result for $\chi$. In
fact, taking into account that at the global minimum
$\partial\Omega/\partial\Delta_{L}=\partial\Omega/\partial\Delta_{R}=0$, and
according to the results in Fig. 3 we have
$\partial\Delta_{L}/\partial\theta=\partial\Delta_{R}/\partial\theta=0$ at
$\theta=0$ (see also Eq. (47)), it is straightforward to prove that (see
appendix B)

$\chi=\frac{\Delta_{L}^{2}}{2G_{D}}\zeta\frac{2-\zeta}{2+\zeta},$ (52)

where $\Delta_{L}$ corresponds to the solution of the gap equation at
$\theta=0$; Eq. (52) stands both at zero and at finite temperature. In Fig. 6
we plot the fourth root of the topological susceptibility versus $\zeta$,
obtained by Eq. (52): we used $\Delta_{L}=50$ MeV and $G_{D}\Lambda^{2}=4.69$:
the agreement between the analytical result and the numerical one is self-
explanatory (we checked the agreement also for other values of the
parameters).

###  III.4 Axion mass and self-coupling

The low energy lagrangian density of the axion field can be written as

$\mathcal{L}_{a}=\frac{1}{2}\partial^{\mu}\partial_{\mu}a-\frac{m_{a}^{2}}{2}a^{2}-\frac{\lambda_{a}}{4!}a^{4},$
(53)

where the axion mass and the quartic coupling are defined in terms of the
potential $V(\theta)=\Omega(\theta)-\Omega(0)$ as

$m_{a}^{2}=\frac{1}{f_{a}^{2}}\left.\frac{d^{2}\Omega}{d\theta^{2}}\right|_{\theta=0},~{}~{}~{}\lambda_{a}=\frac{1}{f_{a}^{4}}\left.\frac{d^{4}V(\theta)}{d\theta^{4}}\right|_{\theta=0},$
(54)

and $\theta=a/f_{a}$. Within our model, we can compute how the axion mass and
coupling behave in the superconductive phase of QCD, as well as study their
response to the phase transition to the normal phase.

From (68) we get an analytical formula for the axion-squared mass in the
color-superconductive phase, which is

$m_{a}^{2}=\frac{\Delta_{L}^{2}}{2G_{D}f_{a}^{2}}|\zeta|\frac{2-|\zeta|}{2+|\zeta|},$
(55)

where $\Delta_{L}$ corresponds to the solution of the gap equation. To our
knowledge, Eq. (55) is a new result in the literature.

Figure 7: Axion self-coupling, $\lambda_{a}f_{a}^{4}$, versus $\zeta$ for
several values of $\Delta_{L}$.

We have not been able to find an analytical expression for $\lambda_{a}$,
hence we computed it numerically. In Fig. 7 we plot the rescaled axion self-
coupling, versus $\zeta$, for several values of $\Delta_{L}$. The parameters
are the same as used in Fig. 6. For the sake of comparison, we note that
studies based on the NJL model find $\lambda_{a}f_{a}^{4}=-(55.64$ MeV$)^{4}$
at $T=\mu=0$ [59]. We find that $\lambda_{a}$ is negative in a range of
$\zeta$ that partly depends on $\Delta_{L}$, hence on the strength of the
coupling. This sign of $\lambda_{a}$ is in agreement with what was found
within NJL models at zero as well as finite $\mu$ [59, 77]. However, we find
also a range of $\zeta$ where $\lambda_{a}$ is positive. The value
$\bar{\zeta}$ of $\zeta\neq 0$ such that $\lambda_{a}=0$ depends on the
strength of the coupling: however, comparing with the results shown in Fig. 3,
we find that $\bar{\zeta}$ is in agreement with the value at which
$\Delta_{L}$ and $\Delta_{R}$ invert their tendency to change as $a/f_{a}$ is
increased, see Eq. (49).

###  III.5 Finite temperature

Figure 8: Critical temperature, $T_{c}$, as a function of $\zeta$, for
different values of $\Delta_{L}$. Calculations correspond to $a=0$.

In this section, we briefly investigate the axion potential in dense
superconductive quark matter at a small, albeit finite, temperature.
Increasing the temperature we expect a critical value, depending on the
chemical potential, for which a phase transition occurs from the color-
superconducting phase to normal-quark matter. In particular, when the critical
temperature $T_{c}$ is reached, the gaps $\Delta_{L}$ and $\Delta_{R}$ vanish
and a second-order phase transition occurs. We do not push this investigation
too much, because we neglected quark masses in this work, that contribute to
obtain nonzero values of $\chi$ and $\lambda_{a}$ above the critical
temperature: we limit ourselves to compute how $T_{c}$ depends on $\zeta$, as
well as to show the behavior of $\chi$ and $\lambda_{a}$ around $T_{c}$
obtained within our model.

In Fig. 8 we plot $T_{c}$ versus $\zeta$ for several values of $\Delta_{L}$;
$G_{D}$ is varied consequently to $\zeta$, such that along each line the value
of $\Delta_{L}$ at $T=0$ is kept constant. We find that $T_{c}$ increases
linearly as a function of $\zeta$; moreover, $T_{c}$ increases upon increasing
$\Delta_{L}$ as expected. This is in agreement with the analytic result that
can be found within the framework of the HDET,

$T_{c}=\Delta_{L}^{0}\frac{e^{\gamma}(2+\zeta)}{2\pi}\,,$ (56)

where $\Delta_{L}^{0}$ is the value of the gap at vanishing temperature and
$\gamma\approx 0.577$ is the Euler-Mascheroni constant. See Appendix C for a
derivation of (56).

It is also interesting to check the behavior of $\chi^{1/4}$ and $\lambda_{a}$
versus $T$ near the transition to normal quark matter. In particular, we
verified that Eq. 52 is also valid in the finite-temperature case (see
appendix B), and the temperature dependence enters only via $\Delta_{L}$. In
the upper panel of Fig. 9 we show the topological susceptibility,
$\chi^{1/4}$ (upper panel) and axion self-coupling, $\lambda_{a}f_{a}^{4}$
(lower panel), versus $T$ for $\zeta=1$. We used $\mu=400$ MeV and
$\Delta_{L}=25$ MeV. The trend we find is in agreement with Eq. 52, with a
vanishing $\chi$ and $\lambda_{a}$ in the normal phase. However, we remark
that the vanishing of these two quantities in the normal phase is an artifact
of our one-loop approximation to the thermodynamic potential, as well as
neglecting the quark masses: including quark masses would lead to nonzero
axion mass and self-coupling also in the normal phase, similarly to what
happens in the high temperature/low density region of the QCD phase diagram
[59, 77, 8].

Figure 9: Topological susceptibility, $\chi^{1/4}$ (upper panel) and axion
self-coupling, $\lambda_{a}f_{a}^{4}$ (lower panel), versus $T$ for $\zeta=1$.
We used $\mu=400$ MeV, $\Lambda=1000$ MeV and $\Delta_{L}=25$ MeV.

##  IV Conclusions and outlook

We analyzed, for the first time in the literature, the QCD axion potential in
dense, superconductive quark matter at finite quark chemical potential, $\mu$,
and finite temperature, $T$. We assumed that both scalar,
$\Delta_{R}-\Delta_{L}$, and pseudoscalar, $\Delta_{R}+\Delta_{L}$, diquark
condensates can form: the color-flavor structure that we assume is that of the
standard 2SC phase. We used a two-flavor model-based on an interaction that
contains a $U(1)_{A}$-preserving and a $U(1)_{A}$-breaking term: the former
has coupling strength $G_{D}$ and is assumed to be derived from one-gluon
exchange, while the latter arises from the one-instanton exchange, it has
coupling strength $\zeta G_{D}$, and contains the coupling of the QCD axion,
$a$, to the quarks. We treated $G_{D}$ and $\zeta$ as free parameters: in
particular, we tuned $G_{D}$ for a given $\zeta$ in order to reproduce a value
of the superconductive gap without axions; we then switch on the axion field
and study the response of the diquark condensate, as well as of the
thermodynamic potential, to $a$. This allowed us to compute the axion
potential in dense quark matter.

We found that for $a=0$ and $\zeta>0$ the scalar condensate is favored;
however, increasing $a/f_{a}$ in the range $(0,\pi/2)$ results in the drying
out of the global minima of $\Omega$ along the direction of the scalar
condensate and in the formation of new global minima around the direction of
the pseudoscalar. Consequently, for $a/f_{a}$ in the range $(\pi/2,3\pi/2)$
the pseudoscalar condensate replaces the scalar one. Increasing further
$a/f_{a}$ up to $2\pi$ the location of the global minima of $\Omega$ change
again and the system turns to condensate in the scalar channel.

We computed the full topological susceptibility, $\chi$, of the 2SC phase,
finding an analytical formula that connects $\chi$ to the quark condensate;
this relation holds both at zero and at finite temperature. Our formula (52)
explicitly contains information relative to the specific model used in our
work, in particular, the effective coupling $G_{D}$ and $\zeta$. It is likely
that in a more refined framework, in which one uses a momentum-dependent one-
gluon exchange term and an instanton kernel instead of our effective
$4-$fermion interactions, both $G_{D}$ and $\zeta$ will be replaced by
quantities directly related to QCD, namely the QCD coupling, a dressed gluon
mass as well as the instanton size. This interesting improvement will be the
the subject of a near future work.

We then computed the axion mass, $m_{a}$, which is related to $\chi$ by the
relation $\chi=m_{a}^{2}f_{a}^{2}$: as a consequence of our result (52), we
were able to obtain an analytical formula for $m_{a}$ in the 2SC phase of QCD.
Also, in this case, it is our hope that the dependence of $m_{a}$ on the
parameters of our model can be replaced by a dependence on quantities of QCD
by using more refined interactions instead of the $4-$fermion terms we adopted
here. We completed this part of the study by computing the axion self-
coupling, $\lambda_{a}$. Interestingly, we found a range of $\zeta$ where
$\lambda_{a}>0$, differently from what it was found in cases of normal quark
matter, see for example [49, 60, 61, 59, 77, 8]: hence the interaction among
axions becomes repulsive in this range of $\zeta$. This might have some impact
on the pressure of axions trapped in the core of compact stellar objects. It
will be interesting to explore this scenario in detail in the future.

A possible improvement of the present work is the inclusion of the strange
quark, and the opening to the three-flavor color superconductor with massive
quarks. This would potentially make the picture more complicated, because in
this work we neglected the current quark masses, as well as light quarks
chiral condensate: previous studies based on NJL-like interactions show that
this can be a fairly good approximation, as long as the light-quarks sector is
considered. However, this might be no longer true for the strange quark, at
least for the values of $\mu$ which are relevant for compact stars. Along this
line, studying the axions in gapless phases is also worth more investigation:
works in this direction are already ongoing and we plan to report on them
soon. Secondly, it would be interesting to use the results obtained here to
compute the cooling of compact stellar objects via axion emission. It would
also be important to check how the picture we drew in our work changes when
the local interaction kernels we used are replaced by non-local ones, with a
more direct link to QCD, as well as when other interaction channels, in primis
the vector and axial-vector channels are included in the game. Even more, it
is of a certain interest to allow for chiral condensate besides the diquark
one, in order to study, within a single model, the axion potential in
proximity of the phase transition between the chiral and the superconductive
phases: this would allow for the exploration of the properties of the axion in
the whole QCD phase diagram. We leave all these interesting improvements to
future works.

###### Acknowledgements.

M. R. acknowledges John Petrucci for inspiration. D.E.A.C. acknowledges
support from the NCN OPUS Project No. 2018/29/B/ST2/02576. A.G. G. would like
to acknowledge the support received from CONICET (Argentina) under Grants No.
PIP 22-24 11220210100150CO and from ANPCyT (Argentina) under Grant No.
PICT20-01847. This work has been partly funded by the European Union – Next
Generation EU through the research grant number P2022Z4P4B “SOPHYA -
Sustainable Optimised PHYsics Algorithms: fundamental physics to build an
advanced society” under the program PRIN 2022 PNRR of the Italian Ministero
dell’Università e Ricerca (MUR).

##  Appendix A HDET gap equation

In this section, we derive the HDET solution to the gap equation (41). To this
end, we put

$\Theta(a,\zeta)=\frac{A(a,\zeta)}{2+\zeta\cos(a/f_{a})}.$ (57)

Then, Eq. (41) gives

$1=\frac{G_{D}}{2\pi^{2}}\Theta(a,\zeta)\int_{0}^{\Lambda}p^{2}dp~{}\left(\frac{1}{\sqrt{(p-\mu)^{2}+G_{D}^{2}h_{L}^{2}A(a,\zeta)}}+\frac{1}{\sqrt{(p+\mu)^{2}+G_{D}^{2}h_{L}^{2}A(a,\zeta)}}\right).$
(58)

We then adopt the approximations of HDET [85]. Firstly, we note that the the
first integral on the right-hand side of (58) gets its largest contribution
from the momentum space region $p\approx\mu$, namely around the Fermi surface
of the quarks. Moreover, the second integral in the right-hand side is
suppressed at large $\mu$ in comparison with the first one, since it does not
receive the enhancement for $p\approx\mu$. Within the spirit of HDET we can
thus ignore the second integral, and restrict the first integral to a thin
shell around $p=\mu$: we call $\delta$ the width of this shell. In the limit
$\delta\ll\mu$ the volume of momentum space available for pairing is thus
$8\pi\mu^{2}\delta$. Introducing $\xi=p-\mu$, the HDET version of (58)
reads111To be precise, in the HDET one introduces a sum over the direction of
the Fermi velocities, $\bm{v}_{F}$, of the quarks; then, for each
$\bm{v}_{F}$, $\xi=\bm{p}\cdot\bm{v}_{F}-\mu v_{F}$ measures the fluctuation
of the longitudinal momentum around the Fermi surface. Our $\xi$ in (59) is
slightly different from that of HDET, nevertheless, formally the gap equation
obtained within HDET is in agreement with (59), because the quark condensate
is homogeneous and the sum over velocities lead to an overall $1$.

$1=\frac{G_{D}\mu^{2}}{2\pi^{2}}\Theta(a,\zeta)\int_{-\delta}^{+\delta}\frac{d\xi}{\sqrt{\xi^{2}+G_{D}^{2}h_{L}^{2}A(a,\zeta)}}.$
(59)

Integration can be done exactly; in the weak coupling limit $\delta\gg
G_{D}h_{L}$, using the definition (39), we finally get

$\Delta_{L}=\frac{4\delta}{\sqrt{A(a,\zeta)}}\exp\left(-\frac{\pi^{2}}{\mu^{2}G_{D}\Theta(a,\zeta)}\right).$
(60)

##  Appendix B Derivation of Eq. (52)

In this section, we derive Eq. (52) starting from $\Omega$ in (37) and using
the gap equation (41). The first step in the calculation is to notice that all
the results we found are consistent with the conditions
$\partial\Omega/\Delta_{L}=\partial\Omega/\Delta_{R}=0$ at $\theta=0$ and at
the minima of $\Omega$. This can be easily understood since $\Omega$ develops
minima along the directions $\Delta_{L}=\pm\Delta_{R}$, and on these lines, it
is an even function of $\Delta_{L}$ or $\Delta_{R}$. Hence we can write

$\chi=\left.\frac{d^{2}\Omega}{d\theta^{2}}\right|_{\theta=0}=\left.\frac{\partial^{2}\Omega}{\partial\theta^{2}}\right|_{\theta=0}.$
(61)

From Eq. (37) we get

$\chi=2\zeta
G_{D}h_{L}h_{R}-\frac{1}{\pi^{2}}\int_{0}^{\Lambda}p^{2}dp\frac{\partial^{2}}{\partial\theta^{2}}\left(\varepsilon_{3}+\varepsilon_{5}+\mu\rightarrow-\mu\right).$
(62)

Along the line $h_{L}=-h_{R}$, which is the relevant one for small $\theta$,
and taking into account the expressions of the dispersion laws of the quarks,
we thus get

$\chi=8\zeta
G_{D}h_{L}^{2}\left[-\frac{1}{4}+\frac{G_{D}}{2\pi^{2}}\int_{0}^{\Lambda}p^{2}dp\left(\frac{1}{\sqrt{(p-\mu)^{2}+(2+\zeta)^{2}G_{D}^{2}h_{L}^{2}}}+\mu\rightarrow-\mu\right)\right].$
(63)

Now, we notice from the gap equation (41) at $\theta=0$ that

$\frac{1}{2+\zeta}=\frac{G_{D}}{2\pi^{2}}\int_{0}^{\Lambda}p^{2}dp\left(\frac{1}{\sqrt{(p-\mu)^{2}+(2+\zeta)^{2}G_{D}^{2}h_{L}^{2}}}+\mu\rightarrow-\mu\right),$
(64)

where $h_{L}$ stands for the condensate at $\theta=0$. Using (64) in (63) we
get

$\chi=2G_{D}h_{L}^{2}\zeta\frac{2-\zeta}{2+\zeta}.$ (65)

Finally, taking into account $\Delta_{L}=2G_{D}h_{L}$ (see Eq. (39)) we have

$\chi=\frac{\Delta_{L}^{2}}{2G_{D}}\zeta\frac{2-\zeta}{2+\zeta},$ (66)

in agreement with Eq. (52). We remark that $\Delta_{L}$ in (66) denotes the
diquark condensate at $\theta=0$. From the derivation presented here, it is
evident that Eq. (52) stands both at zero and at finite temperature, since the
momentum integrals in Eqs. (63) and (64) are modified in the same fashion at
$T\neq 0$.

As we remarked in the main text, for $\zeta<0$ the role of the scalar and
pseudoscalar condensates invert, since the minima of $\Omega$ develop along
the line $h_{L}=h_{R}$ in this case, and Eq. (52) becomes

$\chi=-\frac{\Delta_{L}^{2}}{2G_{D}}\zeta\frac{2+\zeta}{2-\zeta}.$ (67)

Hence, we can summarize the results (52) and (67) as

$\chi=\frac{\Delta_{L}^{2}}{2G_{D}}|\zeta|\frac{2-|\zeta|}{2+|\zeta|}.$ (68)

##  Appendix C Derivation of Eq. (56)

We consider the HDET gap equation Eq. 59, and evaluate the integral
analytically in the case $a=0$ and $\delta\gg G_{D}h_{L}$ thus obtaining

$1=\frac{G_{D}\mu^{2}}{\pi^{2}}(2+\zeta)\ln\left(\frac{4\delta}{G_{D}h_{L}^{0}(2+\zeta)}\right),$
(69)

where $h_{L}^{0}$ is the value of the gap which satisfies the 0-temperature
gap equation.

In the same framework, we now consider the finite temperature gap equation for
$a=0$, which can be written as

$-1+\frac{G_{D}\mu^{2}}{2\pi^{2}}(2+\zeta)\int_{-\delta}^{+\delta}\frac{d\xi}{\sqrt{\xi^{2}+G_{D}^{2}h_{L}^{2}(2+\zeta)^{2}}}=\frac{G_{D}\mu^{2}}{2\pi^{2}}(2+\zeta)\int_{-\delta}^{+\delta}\frac{d\xi}{\sqrt{\xi^{2}+G_{D}^{2}h_{L}^{2}(2+\zeta)^{2}}\left(1+e^{\sqrt{G_{D}h_{L}(2+\zeta)^{2}+\xi^{2}}/T}\right)}.$
(70)

The integral appearing on the left-hand side of the equation above is
analogous to the one in Eq. 59. Thus, the integration leads to the same
result that appears in the r.h.s of Eq. 69, with the only difference that
$h_{L}^{0}$ is replaced by the finite temperature gap $h_{L}$. Inserting then
Eq. 69 into the l.h.s. of Eq. 70, and noticing that the fast convergence of
the integral in the right-hand side of (70) allows us to extend the
integration to the whole real axis, we obtain

$\ln{\left(\frac{h_{L}^{0}}{h_{L}}\right)}=\int_{-\infty}^{+\infty}\frac{d\xi}{\sqrt{\xi^{2}+G_{D}^{2}h_{L}^{2}(2+\zeta)^{2}}\left(1+e^{\sqrt{G_{D}h_{L}(2+\zeta)^{2}+\xi^{2}}/T}\right)}\equiv
I(u),$ (71)

where we put

$I(u)=\int_{-\infty}^{+\infty}\frac{dx}{\sqrt{x^{2}+u^{2}}\left(1+e^{\sqrt{u^{2}+x^{2}}}\right)},$
(72)

with $x=\xi/T$ and $u=G_{D}h_{L}(2+\zeta)/T$. For $T\rightarrow T_{c}$ the
condensate $h_{l}\to 0$: we can thus limit ourselves to consider the leading
order expansion of Eq. 72 for $u\sim 0$, namely

$I(u)\Big{|}_{u\to 0}\simeq\ln\left(\frac{\pi}{e^{\gamma}u}\right),$ (73)

where $\gamma\approx 0.577$ is the Euler-Mascheroni constant. Using this
result in Eq. 71 we obtain

$\ln{\left(\frac{h_{L}^{0}}{h_{L}}\right)}=\ln{\left(\frac{\pi T_{c}}{\gamma
G_{D}(2+\zeta)h_{L}}\right)},$ (74)

where we set $T=T_{c}$ and discarded all the finite $h_{L}^{2}$ contributions.
Finally, Eq. 74 can be fulfilled only if

$h_{L}^{0}=\frac{\pi T_{c}}{\gamma G_{D}(2+\zeta)},$ (75)

which gives from which we trivially derive

$T_{c}=\Delta_{L}^{0}\frac{e^{\gamma}(2+\zeta)}{2\pi},$ (76)

in agreement with (56).

## References

  * [1] R. D. Peccei and H. R. Quinn, Phys. Rev. D 16 (1977), 1791-1797 doi:10.1103/PhysRevD.16.1791
  * [2] R. D. Peccei and H. R. Quinn, Phys. Rev. Lett. 38, 1440-1443 (1977) doi:10.1103/PhysRevLett.38.1440
  * [3] S. Weinberg, Phys. Rev. Lett. 40, 223-226 (1978) doi:10.1103/PhysRevLett.40.223
  * [4] F. Wilczek, Phys. Rev. Lett. 40, 279-282 (1978) doi:10.1103/PhysRevLett.40.279
  * [5] J. E. Kim and G. Carosi, Rev. Mod. Phys. 82, 557-602 (2010) [erratum: Rev. Mod. Phys. 91, no.4, 049902 (2019)] doi:10.1103/RevModPhys.82.557 [arXiv:0807.3125 [hep-ph]].
  * [6] D. A. Easson, I. Sawicki and A. Vikman, JCAP 11, 021 (2011) doi:10.1088/1475-7516/2011/11/021 [arXiv:1109.1047 [hep-th]].
  * [7] E. Berkowitz, M. I. Buchoff and E. Rinaldi, Phys. Rev. D 92, no.3, 034507 (2015) doi:10.1103/PhysRevD.92.034507 [arXiv:1505.07455 [hep-ph]].
  * [8] G. Grilli di Cortona, E. Hardy, J. Pardo Vega and G. Villadoro, JHEP 01, 034 (2016) doi:10.1007/JHEP01(2016)034 [arXiv:1511.02867 [hep-ph]].
  * [9] L. D. Duffy and K. van Bibber, New J. Phys. 11, 105008 (2009) doi:10.1088/1367-2630/11/10/105008 [arXiv:0904.3346 [hep-ph]].
  * [10] M. S. Turner and F. Wilczek, Phys. Rev. Lett. 66, 5-8 (1991) doi:10.1103/PhysRevLett.66.5
  * [11] L. Visinelli and P. Gondolo, Phys. Rev. D 80, 035024 (2009) doi:10.1103/PhysRevD.80.035024 [arXiv:0903.4377 [astro-ph.CO]].
  * [12] P. Sikivie, Phys. Rev. Lett. 48 (1982), 1156-1159 doi:10.1103/PhysRevLett.48.1156
  * [13] J. E. Kim, Phys. Rev. Lett. 43 (1979), 103 doi:10.1103/PhysRevLett.43.103
  * [14] M. A. Shifman, A. I. Vainshtein and V. I. Zakharov, Nucl. Phys. B 166 (1980), 493-506 doi:10.1016/0550-3213(80)90209-6
  * [15] M. Dine, W. Fischler and M. Srednicki, Phys. Lett. B 104 (1981), 199-202 doi:10.1016/0370-2693(81)90590-6
  * [16] A. Caputo and G. Raffelt, doi:10.22323/1.454.0041 [arXiv:2401.13728 [hep-ph]].
  * [17] R. J. Crewther, P. Di Vecchia, G. Veneziano, and E. Witten. Phys. Lett. B, 88, 123 (1979). doi:10.1016/0370-2693(79)90128-X
  * [18] C. A. Baker, D. D. Doyle, P. Geltenbort, K. Green, M. G. D. van der Grinten, P. G. Harris, P. Iaydjiev, S. N. Ivanov, D. J. R. May and J. M. Pendlebury, et al. Phys. Rev. Lett. 97, 131801 (2006) doi:10.1103/PhysRevLett.97.131801 [arXiv:hep-ex/0602020 [hep-ex]].
  * [19] W. C. Griffith, M. D. Swallows, T. H. Loftus, M. V. Romalis, B. R. Heckel and E. N. Fortson, Phys. Rev. Lett. 102, 101601 (2009) doi:10.1103/PhysRevLett.102.101601 [arXiv:0901.2328 [physics.atom-ph]].
  * [20] R. H. Parker, M. R. Dietrich, M. R. Kalita, N. D. Lemke, K. G. Bailey, M. Bishof, J. P. Greene, R. J. Holt, W. Korsch and Z. T. Lu, et al. Phys. Rev. Lett. 114, no.23, 233002 (2015) doi:10.1103/PhysRevLett.114.233002 [arXiv:1504.07477 [nucl-ex]].
  * [21] B. Graner, Y. Chen, E. G. LindahL and B. R. Heckel, Phys. Rev. Lett. 116, no.16, 161601 (2016) [erratum: Phys. Rev. Lett. 119, no.11, 119901 (2017)] doi:10.1103/PhysRevLett.116.161601 [arXiv:1601.04339 [physics.atom-ph]].
  * [22] N. Yamanaka, T. Yamada, E. Hiyama and Y. Funaki, Phys. Rev. C 95, no.6, 065503 (2017) doi:10.1103/PhysRevC.95.065503 [arXiv:1603.03136 [nucl-th]].
  * [23] F. K. Guo, R. Horsley, U. G. Meissner, Y. Nakamura, H. Perlt, P. E. L. Rakow, G. Schierholz, A. Schiller and J. M. Zanotti, Phys. Rev. Lett. 115, no.6, 062001 (2015) doi:10.1103/PhysRevLett.115.062001 [arXiv:1502.02295 [hep-lat]].
  * [24] T. Bhattacharya, V. Cirigliano, R. Gupta, H. W. Lin and B. Yoon, Phys. Rev. Lett. 115, no.21, 212002 (2015) doi:10.1103/PhysRevLett.115.212002 [arXiv:1506.04196 [hep-lat]].
  * [25] M. Colpi, S. L. Shapiro and I. Wasserman, Phys. Rev. Lett. 57, 2485-2488 (1986) doi:10.1103/PhysRevLett.57.2485
  * [26] I. I. Tkachev, Phys. Lett. B 261, 289-293 (1991) doi:10.1016/0370-2693(91)90330-S
  * [27] E. W. Kolb and I. I. Tkachev, Phys. Rev. Lett. 71, 3051-3054 (1993) doi:10.1103/PhysRevLett.71.3051 [arXiv:hep-ph/9303313 [hep-ph]].
  * [28] P. H. Chavanis, Phys. Rev. D 84, 043531 (2011) doi:10.1103/PhysRevD.84.043531 [arXiv:1103.2050 [astro-ph.CO]].
  * [29] F. S. Guzman and L. A. Urena-Lopez, Astrophys. J. 645, 814-819 (2006) doi:10.1086/504508 [arXiv:astro-ph/0603613 [astro-ph]].
  * [30] J. Barranco and A. Bernal, Phys. Rev. D 83, 043525 (2011) doi:10.1103/PhysRevD.83.043525 [arXiv:1001.1769 [astro-ph.CO]].
  * [31] E. Braaten, A. Mohapatra and H. Zhang, Phys. Rev. Lett. 117, no.12, 121801 (2016) doi:10.1103/PhysRevLett.117.121801 [arXiv:1512.00108 [hep-ph]].
  * [32] S. Davidson and T. Schwetz, Phys. Rev. D 93, no.12, 123509 (2016) doi:10.1103/PhysRevD.93.123509 [arXiv:1603.04249 [astro-ph.CO]].
  * [33] J. Eby, M. Leembruggen, P. Suranyi and L. C. R. Wijewardhana, JHEP 12, 066 (2016) doi:10.1007/JHEP12(2016)066 [arXiv:1608.06911 [astro-ph.CO]].
  * [34] T. Helfer, D. Marsh, K. Clough, M. Fairbairn, E. Lim, and R. Becerril, J. Cosmol. Astropart. Phys. 03 (2017) 055 doi:10.1088/1475-7516/2017/03/055 [arXiv:1609.04724 [astro-ph.CO]].
  * [35] D. G. Levkov, A. G. Panin and I. I. Tkachev, Phys. Rev. Lett. 118, no.1, 011301 (2017) doi:10.1103/PhysRevLett.118.011301 [arXiv:1609.03611 [astro-ph.CO]].
  * [36] J. Eby, M. Leembruggen, P. Suranyi and L. C. R. Wijewardhana, JHEP 06, 014 (2017) doi:10.1007/JHEP06(2017)014 [arXiv:1702.05504 [hep-ph]].
  * [37] L. Visinelli, S. Baum, J. Redondo, K. Freese and F. Wilczek, Phys. Lett. B 777, 64-72 (2018) doi:10.1016/j.physletb.2017.12.010 [arXiv:1710.08910 [astro-ph.CO]].
  * [38] P. H. Chavanis, Phys. Rev. D 94, no.8, 083007 (2016) doi:10.1103/PhysRevD.94.083007 [arXiv:1604.05904 [astro-ph.CO]].
  * [39] E. Cotner, Phys. Rev. D 94, no.6, 063503 (2016) doi:10.1103/PhysRevD.94.063503 [arXiv:1608.00547 [astro-ph.CO]].
  * [40] Y. Bai, V. Barger and J. Berger, JHEP 12, 127 (2016) doi:10.1007/JHEP12(2016)127 [arXiv:1612.00438 [hep-ph]].
  * [41] P. Sikivie and Q. Yang, Phys. Rev. Lett. 103, 111301 (2009) doi:10.1103/PhysRevLett.103.111301 [arXiv:0901.1106 [hep-ph]].
  * [42] P. H. Chavanis, Phys. Rev. D 98, no.2, 023009 (2018) doi:10.1103/PhysRevD.98.023009 [arXiv:1710.06268 [gr-qc]].
  * [43] J. W. Wang, X. J. Bi, R. M. Yao and P. F. Yin, Phys. Rev. D 103, no.11, 115021 (2021) doi:10.1103/PhysRevD.103.115021 [arXiv:2101.02585 [hep-ph]].
  * [44] F. Calore, P. Carenza, C. Eckner, T. Fischer, M. Giannotti, J. Jaeckel, K. Kotake, T. Kuroda, A. Mirizzi and F. Sivo, Phys. Rev. D 105, no.6, 063028 (2022) doi:10.1103/PhysRevD.105.063028 [arXiv:2110.03679 [astro-ph.HE]].
  * [45] A. Sedrakian, Phys. Rev. D 93, no.6, 065044 (2016) doi:10.1103/PhysRevD.93.065044 [arXiv:1512.07828 [astro-ph.HE]].
  * [46] L. B. Leinson, JCAP 08, 031 (2014) doi:10.1088/1475-7516/2014/08/031 [arXiv:1405.6873 [hep-ph]].
  * [47] A. Sedrakian, Phys. Rev. D 99, no.4, 043011 (2019) doi:10.1103/PhysRevD.99.043011 [arXiv:1810.00190 [astro-ph.HE]].
  * [48] M. Buschmann, C. Dessert, J. W. Foster, A. J. Long and B. R. Safdi, Phys. Rev. Lett. 128, no.9, 091102 (2022) doi:10.1103/PhysRevLett.128.091102 [arXiv:2111.09892 [hep-ph]].
  * [49] B. S. Lopes, R. L. S. Farias, V. Dexheimer, A. Bandyopadhyay and R. O. Ramos, Phys. Rev. D 106, no.12, L121301 (2022) doi:10.1103/PhysRevD.106.L121301 [arXiv:2206.01631 [hep-ph]].
  * [50] C. H. Lenzi, M. Dutra, O. Lourenço, L. L. Lopes and D. P. Menezes, Eur. Phys. J. C 83, no.3, 266 (2023) doi:10.1140/epjc/s10052-023-11416-y [arXiv:2212.12615 [hep-ph]].
  * [51] G. Lucente, P. Carenza, T. Fischer, M. Giannotti and A. Mirizzi, JCAP 12, 008 (2020) doi:10.1088/1475-7516/2020/12/008 [arXiv:2008.04918 [hep-ph]].
  * [52] T. Fischer, P. Carenza, B. Fore, M. Giannotti, A. Mirizzi and S. Reddy, Phys. Rev. D 104, no.10, 103012 (2021) doi:10.1103/PhysRevD.104.103012 [arXiv:2108.13726 [hep-ph]].
  * [53] Y. Nambu and G. Jona-Lasinio, Phys. Rev. 122, 345-358 (1961) doi:10.1103/PhysRev.122.345
  * [54] M. Buballa, Phys. Rept. 407, 205-376 (2005) doi:10.1016/j.physrep.2004.11.004 [arXiv:hep-ph/0402234 [hep-ph]].
  * [55] S. P. Klevansky, Rev. Mod. Phys. 64, 649-708 (1992) doi:10.1103/RevModPhys.64.649
  * [56] T. Hatsuda and T. Kunihiro, Phys. Rept. 247, 221-367 (1994) doi:10.1016/0370-1573(94)90022-1 [arXiv:hep-ph/9401310 [hep-ph]].
  * [57] S. B. Ruester, V. Werth, M. Buballa, I. A. Shovkovy and D. H. Rischke, Phys. Rev. D 72, 034004 (2005) doi:10.1103/PhysRevD.72.034004 [arXiv:hep-ph/0503184 [hep-ph]].
  * [58] D. Blaschke, S. Fredriksson, H. Grigorian, A. M. Oztas and F. Sandin, Phys. Rev. D 72, 065020 (2005) doi:10.1103/PhysRevD.72.065020 [arXiv:hep-ph/0503194 [hep-ph]].
  * [59] Z. Y. Lu and M. Ruggieri, Phys. Rev. D 100, no.1, 014013 (2019) doi:10.1103/PhysRevD.100.014013 [arXiv:1811.05102 [hep-ph]].
  * [60] A. Bandyopadhyay, R. L. S. Farias, B. S. Lopes and R. O. Ramos, Phys. Rev. D 100, no.7, 076021 (2019) doi:10.1103/PhysRevD.100.076021 [arXiv:1906.09250 [hep-ph]].
  * [61] A. Das, H. Mishra and R. K. Mohapatra, Phys. Rev. D 103, no.7, 074003 (2021) doi:10.1103/PhysRevD.103.074003 [arXiv:2006.15727 [hep-ph]].
  * [62] R. K. Mohapatra, A. Abhishek, A. Das and H. Mishra, Springer Proc. Phys. 277, 455-458 (2022) doi:10.1007/978-981-19-2354-8_83
  * [63] Y. Nambu and G. Jona-Lasinio, Phys. Rev. 124, 246-254 (1961) doi:10.1103/PhysRev.124.246
  * [64] M. G. Alford, A. Schmitt, K. Rajagopal and T. Schäfer, Rev. Mod. Phys. 80, 1455-1515 (2008) doi:10.1103/RevModPhys.80.1455 [arXiv:0709.4635 [hep-ph]].
  * [65] I. A. Shovkovy, Found. Phys. 35, 1309-1358 (2005) doi:10.1007/s10701-005-6440-x [arXiv:nucl-th/0410091 [nucl-th]].
  * [66] R. Rapp, T. Schäfer, E. V. Shuryak and M. Velkovsky, Annals Phys. 280, 35-99 (2000) doi:10.1006/aphy.1999.5991 [arXiv:hep-ph/9904353 [hep-ph]].
  * [67] R. Rapp, T. Schäfer, E. V. Shuryak and M. Velkovsky, Phys. Rev. Lett. 81, 53-56 (1998) doi:10.1103/PhysRevLett.81.53 [arXiv:hep-ph/9711396 [hep-ph]].
  * [68] R. Brower, S. Chandrasekharan, J. W. Negele and U. J. Wiese, Phys. Lett. B 560, 64-74 (2003) doi:10.1016/S0370-2693(03)00369-1 [arXiv:hep-lat/0302005 [hep-lat]].
  * [69] Y. Y. Mao et al. [TWQCD], Phys. Rev. D 80, 034502 (2009) doi:10.1103/PhysRevD.80.034502 [arXiv:0903.2146 [hep-lat]].
  * [70] S. Aoki and H. Fukaya, Phys. Rev. D 81, 034022 (2010) doi:10.1103/PhysRevD.81.034022 [arXiv:0906.4852 [hep-lat]].
  * [71] V. Bernard, S. Descotes-Genon and G. Toucas, JHEP 12, 080 (2012) doi:10.1007/JHEP12(2012)080 [arXiv:1209.4367 [hep-lat]].
  * [72] V. Bernard, S. Descotes-Genon and G. Toucas, JHEP 06, 051 (2012) doi:10.1007/JHEP06(2012)051 [arXiv:1203.0508 [hep-ph]].
  * [73] M. A. Metlitski and A. R. Zhitnitsky, Phys. Lett. B 633, 721-728 (2006) doi:10.1016/j.physletb.2006.01.001 [arXiv:hep-ph/0510162 [hep-ph]].
  * [74] S. Borsanyi, Z. Fodor, J. Guenther, K. H. Kampert, S. D. Katz, T. Kawanai, T. G. Kovacs, S. W. Mages, A. Pasztor and F. Pittler, et al. Nature 539, no.7627, 69-71 (2016) doi:10.1038/nature20115 [arXiv:1606.07494 [hep-lat]].
  * [75] S. Aoki et al. [JLQCD], EPJ Web Conf. 175, 04008 (2018) doi:10.1051/epjconf/201817504008 [arXiv:1712.05541 [hep-lat]].
  * [76] C. Bonati, M. D’Elia, M. Mariti, G. Martinelli, M. Mesiti, F. Negro, F. Sanfilippo and G. Villadoro, JHEP 03, 155 (2016) doi:10.1007/JHEP03(2016)155 [arXiv:1512.06746 [hep-lat]].
  * [77] B. Zhang, D. E. A. Castillo, A. G. Grunfeld and M. Ruggieri, Phys. Rev. D 108, no.5, 054010 (2023) doi:10.1103/PhysRevD.108.054010 [arXiv:2304.10240 [hep-ph]].
  * [78] G. Gabadadze and M. A. Shifman, Phys. Rev. D 62, 114003 (2000) doi:10.1103/PhysRevD.62.114003 [arXiv:hep-ph/0007345 [hep-ph]].
  * [79] F. Takahashi, W. Yin and A. H. Guth, Phys. Rev. D 98, no.1, 015042 (2018) doi:10.1103/PhysRevD.98.015042 [arXiv:1805.08763 [hep-ph]].
  * [80] R. L. Workman et al. [Particle Data Group], PTEP 2022 (2022), 083C01 doi:10.1093/ptep/ptac097
  * [81] Y. Nagashima, Beyond the standard model of elementary particle physics, Wiley-VCH, 2014, ISBN 978-3-527-41177-1, 978-3-527-66505-1
  * [82] M. Shifman, Advanced Topics in Quantum Field Theory, Cambridge University Press, 2022, ISBN 978-1-108-88591-1, 978-1-108-84042-2 doi:10.1017/9781108885911
  * [83] G. ’t Hooft, Phys. Rev. D 14, 3432-3450 (1976) [erratum: Phys. Rev. D 18, 2199 (1978)] doi:10.1103/PhysRevD.14.3432
  * [84] G. ’t Hooft, Phys. Rept. 142, 357-387 (1986) doi:10.1016/0370-1573(86)90117-1
  * [85] G. Nardulli, Riv. Nuovo Cim. 25N3, 1-80 (2002) [arXiv:hep-ph/0202037 [hep-ph]].
  * [86] J. Bardeen, L. N. Cooper and J. R. Schrieffer, Phys. Rev. 106, 162 (1957) doi:10.1103/PhysRev.106.162
  * [87] J. Bardeen, L. N. Cooper and J. R. Schrieffer, Phys. Rev. 108, 1175-1204 (1957) doi:10.1103/PhysRev.108.1175
  * [88] R. Gatto and M. Ruggieri, Phys. Rev. D 85, 054013 (2012) doi:10.1103/PhysRevD.85.054013 [arXiv:1110.4904 [hep-ph]].
  * [89] G. Veneziano, Nucl. Phys. B 159, 213-224 (1979) doi:10.1016/0550-3213(79)90332-8
  * [90] P. Di Vecchia and G. Veneziano, Nucl. Phys. B 171, 253-272 (1980) doi:10.1016/0550-3213(80)90370-3
  * [91] R. J. Crewther, Phys. Lett. B 70, 349-354 (1977) doi:10.1016/0370-2693(77)90675-X
  * [92] H. Leutwyler and A. V. Smilga, Phys. Rev. D 46, 5607-5632 (1992) doi:10.1103/PhysRevD.46.5607