The QCD axion at finite density

We show how the properties of the QCD axion change in systems at finite baryonic density, such as neutron stars. At nuclear saturation densities, where corrections can be reliably computed, we find a mild reduction of the axion mass and up to an order of magnitude enhancement in the model-independent axion coupling to neutrons. At moderately higher densities, if realized, meson (kaon) condensation can trigger axion condensation. We also study the axion potential at asymptotically large densities, where the color-superconducting phase of QCD potentially leads to axion condensation, and the mass of the axion is generically several orders of magnitude smaller than in vacuum due to the suppressed instantons. Several phenomenological consequences of the axion being sourced by neutron stars are discussed, such as its contribution to their total mass, the presence of an axionic brane, or axion-photon conversion in the magnetosphere.


JHEP07(2020)221 1 Introduction
The QCD axion is one of the best motivated particles for physics beyond the Standard Model (SM). The Nambu-Goldstone boson [1,2] of a spontaneously broken U(1) PQ Peccei-Quinn symmetry anomalous under QCD [3], the axion makes the effective QCD θ angle unphysical and it allows QCD dynamics [4] to solve its strong-CP problem, thus explaining the absence of CP violation in the interactions of hadrons, in particular the electric dipole moment of the neutron [5,6]. In addition, the axion constitutes a viable and attractive candidate for dark matter [7][8][9].
The phenomenology of the axion is mainly controlled by the axion decay constant f a . The defining coupling of the axion to the QCD topological charge determines its mass m a ≈ m π f π /f a ≈ 5 µeV(10 12 GeV/f a ) 1 as well as its model-independent couplings to SM fields. Depending on the details of the axion UV completion, extra model-dependent derivative couplings can be present, which scale as well as 1/f a . These properties, and its dependence on f a , endow the axion with a very interesting and diverse phenomenology, which is presently being actively investigated and for which there is a vigorous experimental effort, both current and planned, see e.g. [13]. Particularly relevant probes of the axion are astrophysical in nature, leading to a bound f a 10 8 GeV from the cooling rate of the SN1987A supernova (see e.g. [14] and references therein, as well as [15] and [16] which have recently reevaluated and even put into question the feasibility of such a bound) and f a 10 17 GeV from black hole superradiance [17][18][19].
In this paper we contribute to the global effort in the search for the QCD axion by presenting a first comprehensive study of how its properties change with chemical potential (and negligible temperature), of significant relevance for astrophysical systems, in particular those with large baryonic densities such as (proto-)neutron stars.
Indeed, neutron stars (NSs) are the densest stars in the universe, with densities that in their core are expected to go well beyond nuclear-saturation density, n 0 ≈ 0.16 nucleons/fm 3 (i.e. ρ 0 ≈ m n n 0 ≈ (190 MeV) 4 ≈ 3 × 10 14 g/cm 3 ), which is also a density encountered in supernovae [20][21][22]. Given that such densities correspond to Fermi momenta not far from the scale of (strong) QCD dynamics, and in particular above the pion mass, it is clear that finite density corrections can have an important impact on the properties of the axion. As we will show, at such densities, where the QCD chiral Lagrangian at finite chemical potential is still reliable, the properties of the axion receive corrections of order n/Λ χ f 2 π ≈ 20% (n/n 0 ), in particular leading to a reduction of the axion mass (taking Λ χ ≈ 700 MeV to be the chiral Lagrangian cutoff). A notable exception is the modelindependent coupling of the axion to neutrons, or just the neutron coupling in the KSVZ model, a.k.a. hadronic axion [23,24]. We find that the accidental cancellation present in vacuum (see e.g. [25]) is spoiled at finite density, potentially resulting in an O(10) enhancement, up to uncertainties of the relevant operators' coefficients and on yet to be computed extra corrections.

JHEP07(2020)221
Since chiral perturbation theory (ChPT) is expected to fail beyond nuclear-saturation densities, and there is no other known reliable method to describe QCD at such high densities as those found in the core of NSs, we are driven to make extrapolations and qualitative statements about the possible fate of the QCD axion. However, under well defined assumptions, in particular that kaon condensation [26] takes place somewhat above n 0 , we show that spontaneous CP violation in the neutral Goldstone sector is triggered, resulting in axion condensation -see [27] for a first discussion of axion condensation in the context of a tuned non-standard axion. For a range of axion masses of order of the inverse radius of the NS, the axion gets sourced, with phenomenological consequences that go from a non-standard contribution to the mass of the star (closely resembling [28,29]) or the presence of an axionic brane [30], to an O(1) enhancement of the electric field in the star's magnetosphere. Other implications of such a sourcing, e.g. axion-mediated longrange forces with a potentially dramatic impact on NS mergers [27], do not seem to be realistic for the QCD axion.
To improve our understanding, at least at the qualitative level, of the physics of the axion at high densities, we also resort to the limit of very high quark chemical potential, µ q Λ χ . In this regime, QCD is expected to reside in a color-flavor-locked (CFL) phase, where gluons are massive and a perturbative expansion is possible [31][32][33][34]. Our analysis of the axion potential in this regime shows that axion condensation is a generic phenomenon (as indirectly anticipated in [35], where a non-zero effective θ angle was found), while the axion mass is orders of magnitude smaller than in vacuum, mainly due to the fact that QCD-instantons are very suppressed, scaling as 1/µ 8 q . The fact that we find the same qualitative behavior for both the edges of QCD perturbative control, i.e. at mild and very high densities, is very suggestive of the fact that a non-zero QCD-axion condensate could indeed develop in NSs.
We would like to note that our analysis of the physics of a light field such as the QCD-axion in systems at finite chemical potential is of general applicability, therefore of relevance for other light (below the QCD confinement scale) axion-like particles coupled to QCD degrees of freedom, e.g. [17,[36][37][38][39]. 2 We find this type of new physics well motivated and very interesting, specially in light of the excellent experimental prospects for improving our understanding of dense astrophysical objects such as NSs, as exemplified by the recent detection of gravitational waves from merger events [41].
The rest of the paper is organized as follows. We review the axion potential in vacuum in section 2 and how the chemical potential is introduced in quantum field theory, with illustrative examples of boson condensation and two-flavor QCD meson condensation, in section 3. Readers familiar with these topics may skip to section 4, where we present our results for the properties of the axion in nuclear matter around nuclear-saturation density. The study of the axion potential in the CFL phase, and our extrapolation to densities expected to be found in NSs cores, is found in section 5. An overview of the potential 2 A notable scenario is that of the realizations of the relaxion [36] which rely on QCD dynamics as the main source of stabilization (back-reaction), see also e.g. [37]. The potentially dramatic effects of baryon chemical potential and the corresponding implications on NSs for these scenarios will be discussed elsewhere [40].

JHEP07(2020)221
effects of an axion condensate in NSs is presented in section 6, deferring for future work a more detailed study of this and other implications on NSs and supernovae [30]. We present our conclusions and outlook in section 7. Finally, several appendices are devoted to detailed calculations of the results presented in the main text.

Axion potential in vacuum
The most general QCD and axion effective Lagrangian below the electroweak scale, at leading order in fields and derivatives, is given by (2.1a) with implicit flavor and color indices. J µ PQ = q c 0 qq γ µ γ 5 q is a model-dependent current associated with a spontaneously broken axial U(1) PQ symmetry, made of the SM matter fields q. The Nambu-Goldstone boson (NGB) of the U(1) PQ is the axion field a(x), with decay constant f a defined by its coupling to gluons. The axion coupling to photons is given by g 0 aγγ = e 2 8π 2 fa E N , with E/N the ratio of the electromagnetic (EM) and the color anomalies. In the free and chiral limit the theory is invariant under the symmetry group At low energies, QCD confines and a chiral condensate q R q L develops that breaks spontaneously the global symmetries The low energy degrees of freedom are described by the fluctuations of the condensate, i.e. the NGBs of the broken chiral symmetries 3

JHEP07(2020)221
The explicit breaking of the chiral symmetries by the quark masses can be incorporated in the low-energy theory by promoting the quark mass matrix, M = Diag[m u , m d , m s ], to a spurion with the transformation properties Under a U(1) axial rotation, the θ angle shifts as θ → θ + 2N f α A . The U(1) A is therefore anomalous, explicitly broken by non-perturbative effects associated with incalculable large instantons. Since the shift symmetry of the axion, associated with U(1) PQ , can be used to remove the θ angle from the Lagrangian eq. (2.1), a → a − θf a , the axion can be treated as an actual dynamical spurion for the U(1) A . The non-perturbative nature of the axial anomaly means that the effective Lagrangian for the η , which shifts under U(1) A as η → η + 2N f α A f η , is not calculable. 4 That would be the case for the axion as well, if not for the fact that one can move to a different basis by performing a local chiral transformation of the quarks in eq. (2.1), with Q a an arbitrary matrix in flavor space which, if Tr[Q a ] = 1, eliminates the axion coupling to gluons. In this basis, the Lagrangian above the QCD confinement scale reads with Q e = Diag[2/3, −1/3, −1/3] the flavor-space matrix of electric charges. After such a redefinition of the quark fields, and upon integrating out the heavy η , the axion potential can be related to that of the QCD pions where m π is the neutral pion mass and we neglected O(∆m/m s ) terms, ∆m ≡ 1 2 (m u −m d ). In this N f = 2 approximation, the quark condensates evaluated in the vacuum are 2b = qq 0 ≡ 1 2 ūu +dd 0 , which leads to the Gell-Mann-Oakes-Renner (GOR) relation (2.12)

JHEP07(2020)221
The axion mass can be calculated at leading order by integrating out the chiral NGBs at tree-level, see appendix A. The final result for the axion mass reads where we neglected corrections of order O(m u,d /m s ), since they are numerically of the same order as other next-to-leading order (NLO) corrections (e.g. η mixing) [25]. We finally note that, as it could not be otherwise, the axion mass is independent of the arbitrarily chosen Q a . Such a matrix however can be chosen to simplify the calculation of a given observable. For example, choosing Q a = M −1 Tr[M −1 ] removes all the tree-level mixing between the axion and the neutral mesons, thus simplifying e.g. the calculation of the axion mass.
3 Chemical potential in quantum field theory Introducing a chemical potential in quantum field theory is a generalization of the procedure in statistical mechanics. One defines a new operator corresponding to the thermodynamic Landau free energy (a.k.a. grand thermodynamic potential density) with H the Hamiltonian density, J 0 i the conserved charge density associated with a given global symmetry of the system (i.e. the temporal component of the conserved current), and µ i the corresponding chemical potential. 5 From the path integral representation of the partition function (see e.g. [42]), one arrives at the following prescription: the temporal derivative of each field transforming under the global symmetry in question is shifted by with T R i the generator of the global symmetry in the appropriate representation R. Chemical potential therefore acts as a source for the temporal component of the corresponding conserved current, much like a background gauge field potential. Since it singles out the time direction, the chemical potential breaks the Lorentz symmetry down to its SO(3) subgroup of spatial rotations. Charge conjugation symmetry (C), under which J 0 i → −J 0 i , is also broken, while parity (P) and time-reserval (T) are preserved -CP and CPT are thus broken. If part of a non-abelian group, a chemical potential also breaks the global symmetry by singling out a specific direction in generator space, namely µ i T i , which defines an unbroken U(1) subgroup. 5 We recall that the grand-canonical density matrix is given byρ = exp [−β(H − µiQi)], with β = 1/T (T is the temperature), H the Hamiltonian, and Qi the conserved charge. The partition function is then Z(V, T, µi) = Trρ, where V is the volume, eventually taken to infinity. The thermodynamic potential density is Ω(T, µ) = −(T /V) ln Z = ρ − µini = −p, with ρ the energy density, ni the number density, and p the pressure. The grand-canonical average of an operator O is then O T,µ i = Tr[Oρ]/Z (with a slight abuse of notation, when clear we will denote ensemble averages simply by O ). Then ni = J 0 i = −(∂Ω/∂µ)T , while the entropy density is given by s = −(∂Ω/∂T )µ.

U(1) toy model
A simple toy model that illustrates the main effect of the chemical potential is a complex scalar theory with a global U(1) symmetry [43,44]. After using the prescription of eq. (3.2), one finds the following Lagrangian For m 2 > µ 2 , the field expectation value is trivial, φ = 0, and respects the global U(1) symmetry. The two propagating degrees of freedom have different dispersion relations The appearance of the chemical potential breaks C symmetry, which appears as a φ ↔ φ * exchange symmetry in the µ = 0 theory -therefore µ can be treated as a spurion transforming as µ → −µ.
Above the threshold |µ| > m, the global U(1) is spontaneously broken by the expectation value and the theory describes a Bose-Einstein condensate (BEC) phase. In contrast to the ideal (λ = 0) ultra-relativisitic Bose gas [43], in the interacting theory (with λ > 0) the chemical potential can be larger than m [44], without leading to any inconsistencies. Note, that our fundamental potential being the Landau free energy Ω, the fixed thermodynamical parameter is µ, which sets the effective energies of the particles in the system due to a coupling to the "particle bath". This allows the flow of particles in and out of the system, implying that the charge density, n φ − n φ * , is a derived quantity set by µ. 6 As we show below, for |µ| > m, the T = 0 system contains non vanishing charge density in the form of the BEC. One can interpret this appearance of charge as particles from the "particle bath" being inserted in the ground state of the system.
In the BEC phase, the following parameterization is useful The classical potential is minimized for v 2 = µ 2 −m 2 λ , and we find the following Lagrangian . The charge density in the condensed phase is non-vanishing in the limit of zero temperature β ≡ 1/T → ∞ and infinite volume V → ∞, JHEP07(2020)221 where we used the classical ( → 0) result for the generating functional ln Z = −βVV (µ) for a homogeneous classical configuration σ . By diagonalizing the quadratic field operators in momentum space one finds the dispersion relations for the two propagating degrees of freedom which at zero momentum are As expected, there is one massless excitation, corresponding to the NGB of the spontaneously broken U(1), and one massive excitation, the radial (or Higgs) mode.

Meson condensation
We review now the importance of a chemical potential in the context of meson condensation in QCD, in particular for the case of two flavors [45][46][47][48], and discuss for the first time its effects on the axion potential. This is a simplified version of the more complicated, but plausibly more realistic, scenario of kaon condensation (N f = 3), to be discussed in section 4.2.
For N f = 2, the chiral condensate breaking SU(2) L × SU(2) R × U(1) B × U(1) A spontaneously to SU(2) × U(1) B can be parameterized, in full generality, as q R q L ≡ q R q L 0 e iα Σ 0 , Σ 0 = cos θ 1 2 + i sin θn · σ ,n = (sin ψ cos χ, sin ψ sin χ, cos ψ) , (3.10) where −π/2 ≤ θ < π/2, 7 and with the sigma field transforming as In eq. (3.10), we used the fact that a field transforming as a bi-fundamental under SU(2) L × SU(2) R can be written as a radial mode, here frozen to some constant value q R q L 0 , times a 2-by-2 unitary matrix Σ 0 , which parameterizes the orientation of the ensemble average in the presence of finite µ, which we call the orientation of the expectation value here. 8 The phase factor e iα is identified with the direction in field space associated with the anomalous axial U(1). A potential for α is generated by non-perturbative effects, whose minimum is at α = 0, which we take from this point on. The angles defined in eq. (3.10) can be related to expectation values of the usual pion fields (at vanishing chemical potential) The shift θ → θ + π can be compensated by shifting α → α + π. 8 SU(2)L and SU(2)R are generated by T a L = 1 2 σ a and T a R = 1 2 σ a , respectively, where as usual it should be understood that the L and R operators act on different indices and therefore commute.

JHEP07(2020)221
where we defined π i π i ≡ Π . In Dirac notation where we denoted q R q L 0 = q L q R 0 ≡ qq 0 /2. Therefore CP is broken in the ground state if Σ 0 = Σ † 0 , that is if θ = 0. We wish to study this system at a non-vanishing chemical potential for isospin and we shall neglect for the remainder of this section isospin breaking due to the quark masses and electromagnetic interactions, making the choice in eq. (3.14) completely generic. Such a chemical potential is associated with the σ 3 rotation of the vector SU(2) L+R subgroup. Therefore, according to eqs. (3.2) and (3.11), we promote the temporal derivative of Σ 0 to Note that changingμ →μ+ 1 6 1 2 in eq. (3.14) has no effect on eq. (3.15) and on the following derivation, therefore in this context the isospin chemical potential can be equivalently associated with the chemical potential for electric charge. 9 The resulting potential for the pions and the axion, the latter entering via the quark mass matrix, M = m1 2 (m u = m d ), as in eq. (2.9) (with Q a = 1 2 /2), is given by at leading order in m/Λ χ and µ/Λ χ , Λ χ being the cutoff of the chiral Lagrangian. We note that the first term arises from the usual kinetic term, 1 4 f 2 π Tr[∂ µ Σ 0 ∂ µ Σ † 0 ], after the replacement (3.15). Using eq. (3.10) we find where m π here is the neutral pion mass in vacuum, i.e. m 2 π f 2 π = −2m qq 0 . We see that the isospin chemical potential introduces an additional source of explicit symmetry breaking -while leaving unbroken the U(1) L+R symmetry defined by the generator in eq. (3.14),μ explicitly breaks the shift symmetries associated with the would-be NGBs charged under U(1) L+R , i.e. the charged pions. Indeed, as discussed above, U(1) L+R is equivalent to the electric charge. Consequently, the first term in eq. (3.17) is proportional to the expectation value of the charged pions, sin 2 θ sin 2 ψ ∝ π + π − . Sinceμ commutes with the U(1) L−R associated with the neutral pion, the neutral NGBs are unaffected by the chemical potential and the potential eq. (3.17) is minimized at π 3 = 0 (ψ = π/2) and a = 0 as in the µ = 0 vacuum.

JHEP07(2020)221
The minimum of the potential for any value of µ is then found at For |µ| < m π , the ground state is the trivial one, Σ 0 = 1, thus its orientation is the same as for µ = 0. For |µ| > m π , pion condensation takes place and the orientation of the expectation value is no longer trivial. We note that in this case χ constitutes a flat direction which, as we confirm later, corresponds to a NGB from the spontaneous breaking of electric charge, U(1) L+R . Setting, without loss of generality, χ = 0, we can write the QCD orientation for |µ| > m π as Σ 0 = cos θ i sin θ i sin θ cos θ . (3.19) At this point we recall that since θ = 0, CP is broken by the expectation value, a result of a sufficiently large explicit breaking of CP by the chemical potential in the charged pion sector. Instead, CP-invariance in the neutral sector is preserved by the charge chemical potential, which leaves the expectation values in that sector untouched. We see now that only if π 3 = 0 (ψ = π/2) could the axion condense, which requires additionally explicit breaking of isospin, i.e. m u = m d .
Having established the Goldstone boson expectation values at finite-density, let us turn our attention to their fluctuations. Since these are associated with the SU(2) L × SU(2) R generators broken by Σ 0 , we define the following rotated generators where ξ 0 ≡ √ Σ 0 . The broken and unbroken generators are then given by respectively. The fluctuations around the Σ 0 ground state can be parameterized as where, abusing notation, we have written the (pseudo-)NGBs as π a , like the standard θ = 0 pions. 10 The dispersion relations for the neutral degrees of freedom, π 0 and the axion, are the same as for vanishing chemical potential. Their masses can be obtained from eq. (3.16) (with the substitution of Σ 0 by Σ), (m 2 π 0 ) θ = m 2 π / cos θ , (m 2 a ) θ = (m 2 a ) 0 cos θ , (3.23) 10 We note that, given eq. (3.20) and T a L = 1 2 σ a , T a R = 1 2 σ a , it follows that ξL = ξ0 exp iπ a σ a 2fπ ξ † 0 and The masses are normalized to their respective µ = 0 value. The charged π + and π − modes (orange and blue curves respectively) evolve similarly as the φ and φ * modes in the U(1) toy model: a linear split in masses in the uncondensed phase, continuously transitioning to a massless Goldstone mode and a massive radial mode in the condensed phase. The masses of the neutral modes, π 0 and a (green and red curve respectively) are unaffected by the chemical potential in the uncondensed phase. In the condensed phase, m π0 increase linearly with µ, while the axion becomes lighter as µ increases.
with (m 2 a ) 0 the mass of the axion in vacuum, eq. (2.13), and where we note that for |µ| > m π , (m 2 π 0 ) θ = µ 2 . The change of the axion mass for θ = 0 simply follows from the fact that, once the mixing with π 3 is eliminated, it has to be proportional to the CP-even combination Tr[Σ 0 + Σ † 0 ] ∝ cos θ. The increase in the neutral pion mass can be understood as a result of its repulsive interaction with the charged pions. The dispersion relation for the charged pions is very similar to the U(1) toy model of section 3.1. In the uncondensed phase |µ| < m π , their dispersion relations are Indeed, for the charged states π ± ≡ 1 √ 2 (π 1 ∓ iπ 2 ) we recognize the same mass splitting we found in eq. (3.4). In the condensed phase |µ| > m π , the remaining U(1) L+R symmetry is spontaneously broken. The effective masses of the charged pions are As in the U(1) toy model, the condensed phase contains one massless Goldstone mode and one massive radial mode. In this phase, the system has a non-vanishing charge density The effective masses of the pions and the axion are plotted in figure 1.

JHEP07(2020)221 4 Nuclear phase
In this section we study how the properties of the axion, mainly its potential and coupling to nucleons, change in systems at finite baryon density, n. In particular, our focus here is on densities around nuclear saturation, n ∼ n 0 , where a description of QCD in terms of hadrons is still meaningful. For the axion potential, we identify two main effects: 1) the change in the size and, to some degree, flavor orientation of the quark condensates, as "measured" by the mass of the pions (section 4.1), and 2) kaon condensation (section 4.2), similar to meson condensation, introduced in section 3.2. Both of these effects can be taken into account by a generalization of the axion potential in vacuum, eq. (2.10), to with M a encoding the dependence on the axion as in eq. (2.9). Σ 0 = ξ 2 0 parametrizes the orientation of the QCD ground state that spontaneously breaks SU(3) L × SU(3) R to SU (3) and therefore encodes the effects of kaon condensation. In vacuum, we have ξ 0 = 1 3 and the unbroken subgroup is the usual SU(3) L+R , while in the kaon-condensed phase, we have ξ 0 = ξ 0 (θ), with θ controlling the size of the kaon condensate which, as explained below, ultimately depends on the baryon density. ξ L,R are the Goldstone matrices, given by a generalization to SU(3) L × SU(3) R of those in eq. (3.22). Finally, the quark condensate qq n at finite density becomes a matrix in flavor space, 3) The detailed derivation of eq. (4.1) is given in appendix B. The result can also be understood in terms of symmetries:M a is a spurion that has been dressed by the Goldstones and projected into the SU(3) L+R subgroup. Therefore, it transforms asM a → VM a V † , where V is an SU(3) L+R transformation. qq n transforms in the same way, since it is the result of a non-vanishing expectation value of the temporal component of the baryonic current, n = J 0 B . 11 Concerning the couplings of the axion to nucleons, the main effect we consider can be traced to a change at finite baryonic density of the nuclear matrix elements p|qγ µ γ 5 q|p (with p the proton, q = u, d) as "measured" by the axial pion-nucleon coupling (section 4.3). These set the size of the couplings of protons and neutrons to the axion, as they follow from either its model-dependent UV couplings to light quarks, or from axion-pion mixing. The latter also changes at finite density although, as we explain below, the effect is within known uncertainties.

JHEP07(2020)221
Before going into the details, several general comments about our treatment of the nuclear medium are in order. To describe the state of the system, we will work directly in terms of baryon densities, n p and n n considering only protons and neutrons, respectively. In practice, our independent parameters are the total baryonic density, n = n p +n n , and the proton fraction, n p /n. This will be more convenient than introducing the corresponding chemical potentials, because our analysis is limited to linear order in n, i.e. we work in the mean-field or Hartree approximation, where e.g. n p ≈ pγ 0 p T,µ i (and in fact n p ≈ pp T,µ i in the non-relativistic limit) -higher-order corrections generically being beyond perturbative control when relevant. Besides, the relative fraction of protons and neutrons is, as shown below, relevant only in our discussion of kaon condensation. There, the chemical potential for electric charge, µ, will also be required to properly describe the system, along with the condition of charge neutrality.

Quark condensates
We first discuss how the quark condensate changes at finite baryonic density, since this is the most robust effect from the point of view of perturbative control. We derive the implications for the axion mass, which were first noted in [27]. The change with density of the quark condensates can be calculated utilizing the Hellmann-Feynman theorem [49] ∆E(n) is the energy shift of the QCD ground state due to the finite density background, such that ∆E(0) = 0. It can be decomposed as where the first term represents the energy shift due to the presence of a non-interacting Fermi gas, while the second term encodes the energy shift due to nuclear interactions. Neglecting these interactions as well as relativistic corrections, we have ∆E = x=n,p,... m x n x , and we arrive at the so-called linear approximation for the in-medium condensate The derivatives ∂m x /∂m q describe the shift in the nucleon mass due to the non-vanishing quark masses. For two nucleons {n, p} and three quarks {u, d, s}, one naively counts six independent shifts. However, due to the {p, u} ↔ {n, d} exchange symmetry, only three shifts are independent. Working in the isospin basis for the quark masses,m ≡ 1 2 (m u + m d ) and ∆m ≡ 1 2 (m u − m d ), the following sigma terms are identified and defined

JHEP07(2020)221
such that with M B the baryon mass in the chiral limit, m q → 0. We note that the sigma terms can be expressed in terms of the parameters of the N f = 3 chiral Lagrangian for baryons, see eq. (B.38) in appendix B. The σ πN and σ s terms have been extracted from pion-nucleon and kaon-nucleon scattering experiments, as well as from lattice simulations by calculating the mass shifts of the nucleons. There are ongoing efforts aimed at the determination of the precise values of these sigma terms. A summary of latest results [50] shows that their current values are scattered over a fairly wide range, with some tension between experimental and lattice results. In this work we use the conservative estimates σ πN = 45 ± 15 MeV and σ s = 30 MeV . The other sigma term is extracted from the p − n nonelectromagnetic mass splitting 2σ πN = (m n − m p ) non-EM = 2 ± 0.3 MeV [51]. Using the GOR relation in eq. (2.12), we rewrite the ratios qq n / qq 0 as Clearly the ss n condensate is only weakly affected by the nucleonic background. Therefore, as in vacuum, its contribution to the axion mass will be subleading, being suppressed by m u,d /m s . Additionally, ūu n ≈ d d n up to the small isospin breaking correction [52], which we neglect. From eq. (4.1) with ξ 0 = 1 3 and after taking care of axion-pion mixing (which we discuss in the context of the axion couplings section 4.3) we reproduce the axion mass at finite density found in [27] (m a ) 2 n = where m π here is the neutral pion mass in vacuum, eq. (2.12). In this regard, we note that at the linear order in density the same correction as the axion enters the neutral pion mass in medium, i.e. (m π ) 2 n = m 2 π ūu n . This is why for the remainder of this section, we shall only consider n < n c ≡ n 0 /b 1 ≈ 2.8 n 0 (45 MeV/σ πN ), with n c being the critical density in which one naively expects chiral symmetry restoration in the linear approximation.

JHEP07(2020)221
At this point, let us turn our attention to the corrections to the linear, non-relativistic approximation we have considered. This will allow us to estimate the densities up to which our leading result is under perturbative control and can therefore be trusted. First, the energy of a degenerate (zero temperature) ideal Fermi gas receives relativistic corrections. In the fully relativistic limit, the free part of the energy for a fermion x is given by where x , which determines the number density, (4.17) , become important at large densities. When this happens, corrections to the QCD ground state energy eq. (4.5) from nucleon interactions become important as well. These are predominantly due to pion exchange, but also from four-baryon contact interactions. It is clear that the latter become important when n x /Λ χ f 2 π becomes order one. Given eq. (4.17), this is also the place where ChPT is beyond control, k x f ∼ Λ χ , e.g. the pion-exchange contribution to the energy is not predictable. In addition, since the cutoff of ChPT Λ χ is numerically close to m p ≈ m n , relativistic corrections are approximately controlled by the same expansion parameter, where we took k f = k p f ∼ k n f . Ultimately the best way to asses the validity of our linear approximation is to explicitly compute the relevant NLO corrections. The interaction energy E int. has been calculated by summing the so-called Hugenholtz diagrams, which are connected bubble diagrams describing ground-state to ground-state transitions [53]. The resulting higher-order finite density effects on the quark condensates have been obtained in ChPT for symmetric nuclear matter [54,55] and pure neutron matter [56,57]. These authors have indeed found O(1) deviations from the linear approximation for densities somewhat above nuclear saturation. Specifically, nucleon interactions seem to ameliorate the linear decrease of ūu n ≈ d d n in eq. (4.12), such that already at n ≈ 2n 0 , the condensates are only at approximately 60 % of their vacuum value, as opposed to the 15 % predicted by the linear approximation, and in fact start increasing with density [54]. This then implies that a more realistic prediction of the axion mass in dense symmetric nuclear matter is while for larger densities n 2 0 ∼ n c it becomes difficult to trust the results of ChPT. Therefore, the determination of the quark condensates and the axion mass at densities significantly beyond nuclear saturation remains an open and difficult theoretical problem.

JHEP07(2020)221
Importantly, realistic lattice simulations at finite density are currently not feasible due to so-called sign problem. In addition, at such high densities other issues arise (ultimately related to the problem of perturbativity), such as the "hyperon puzzle", which concerns the appearance, or absence, of hyperons, see e.g. [58] and references therein. In the next section we will focus our attention instead on another effect of strangeness, potentially much more relevant for the fate of the axion at finite density.

Kaon condensation
In the previous section we assumed that the vacuum of QCD is trivially oriented and CP-preserving, or equivalently that none of the mesons acquire a non-trivial expectation value. This might, however, not be the case in dense matter. It has been hypothesized [26] that above certain baryonic densities it becomes energetically possible for the strangeness changing process of a neutron splitting into a proton and a scalar K − meson, and vice versa, to take place The reason being the low in-medium kaon mass, which eventually leads to the formation of a K − condensate. This process takes place along with, and even becomes favored over, the usual neutron β-decay, n → p + + e − +ν e , and inverse β-decay, p + + e − → n + ν e , because of the high price of occupying the increasingly energetic Fermi surface of the electrons. Because of this fact, also the processes e − ↔ K − + ν e and e − ↔ µ − +ν µ + ν e , µ − ↔ e − + ν µ +ν e reach β-equilibrium [59]. On the other hand, the formation of a pion (π − ) condensate (n ↔ p + + π − ) seems to be disfavored, as we shall discuss below. Motivated by these arguments, we shall now entertain the possibility of kaon condensation and derive its effects on the axion potential. Several important comments and some caveats are however in order. We consider this scenario because of the thrilling possibility of leading to axion condensation, even though it takes place -if it takes place at all -at densities where a perturbative expansion is questionable, n 2n 0 . Because of the inherent uncertainties at such densities, our conclusions will be qualitative rather than quantitive. Indeed, similar to our discussion at the end of the previous section on the quark condensates and their finite density corrections beyond the linear approximation, kaon condensation cannot be simply described by the leading order terms in ChPT. In particular, nucleon selfinteractions and interactions with pions need to be considered in order to capture the full complexity of this strongly interacting system [60] -for instance, the latter are the reason behind the fact that K − condensation is more likely than π − condensation. Our working assumption is that all the processes above (neglecting the pions) are in equilibrium, which implies a set of equations relating the chemical potentials of the particle species involved, where µ is the chemical potential associated with (positive) electric charge. For convenience, we work directly with muon and electron densities, n µ and n e respectively, both of which are determined by µ as they follow from an ideal Fermi gas. The size of the kaon condensate, θ, is determined, as in the simple example of meson condensation discussed in JHEP07(2020)221 section 3.2, by the minimization of the scalar potential, which of course also determines if the axion condenses or not. Finally, due to the importance of nuclear interactions, the densities of protons and neutron, or equivalently the total baryon density n and the proton fraction n p /n, are not determined by µ. Instead, we enforce the condition of (electric) charge neutrality n EM = 0, and present our results in terms of n and n p /n. An additional important final comment regards the implications of kaon condensation on the NS equation of state (EoS). It has been argued that the inclusion of kaon condensation generically leads to a softer EoS [59,61], which usually cannot sustain a large NS mass. This is in conflict with the most massive NSs observed, with masses around 2M [62,63]. This is the main reason why kaon condensation is currently considered an "exotic" possibility. However, axion effects can in fact harden the EoS [30] and reopen this window. Also, kaon condensation is in fact related to another issue, namely the hyperon puzzle [64]. The appearance of hyperons also tends to soften the EoS, resulting in a similar apparent conflict with the observation of massive NSs. Therefore, although the appearance of strangeness seems to be in tension with observations, we think it would be premature to definitively exclude the possibility of kaon condensation at this point, especially in the presence of new physics.
Let us consider then the possibility that kaon condensation occurs in nuclear matter and qualitatively examine its effects on the axion potential. Once the chemical potential for electric charge µ is introduced, the dispersion relations for the K ± modes are given by with the kaon effective in-medium mass The first term in eq. (4.23) is the usual kaon mass to leading order in m s , with the inclusion of the finite density corrections to the relevant quark condensate, which in the linear approximation are given by where m 2 K = −m s qq 0 /f 2 π , the neutral kaon mass in vacuum, neglecting O(m u,d /m s ) terms. The second term in eq. (4.23) is a mass correction induced by the baryonic background, due to the model-independent s-wave interactions of the baryons with the mesons, arising from the baryon kinetic term, as it follows from the covariant derivative of ChPT,

JHEP07(2020)221
Note that since b 2 b 1 , the effective kaon mass decreases with density, and condensation is expected to occur when 12 Kaon condensation is introduced by allowing the kaon field to take a non-trivial average value, √ 2K ± /f π , or equivalently, in our notation, reorienting the QCD ground state in medium [59], (4.28) The ground state orientation is determined by the static Lagrangian after setting all the fluctuations to zero. Neglecting for the time being the axion, we find a similar potential to eq. (3.17) Minimizing V (θ) leads to the condition which can be used to determined θ = θ(µ, n, n p /n). The requirement of electrical neutrality, n EM = − ∂L/∂µ = 0, leads to − f 2 π µ sin 2 θ + cos θ n p − sin 2 (θ/2) n n − n e (µ) − n µ (µ) = 0 , (4.31) where we included the lepton charge densities, given by  figure 4 we plot the region (blue) in the {n, n p /n} plane where θ = 0, namely where the system is in the kaoncondensed phase. The evolution for given {n, n p /n} can be understood as follows: for a fixed proton fraction n p /n, as n increases the amount of positive charge due to the protons increases as well, and more leptons are required to satisfy the neutrality condition, leading 12 It is illustrative to also consider the pion effective in-medium mass, since it shows that, due to the second term and contrary to the kaon, the charge pion becomes heavier with increasing density, at least for a neutron rich background np/n < 1/2 [65]. The argument against pion condensation becomes even stronger when considering higher order terms in ChPT [52], as we discussed at the end of section 4.1 -these additional corrections even make the pion mass increase with density for n ≈ 2n0, even in symmetric nuclear matter, np/n = 1/2. to an increase in µ. This increase in µ drives the effective mass of the kaon, eq. (4.23), further down (on top of the decrease in ūu at finite density), until eventually the threshold condition for kaon condensation is met, µ = (m K ± ) n , and a further increase in the proton density can be compensated by inserting K − particles in the ground state. Even though we keep them undetermined, let us briefly comment at this point on how the proton fraction could be determined in terms of the total density. Since the interaction energy of nuclear matter also depends on the proton fraction, E int. = n int. (n, n p ), one could enforce that the total energy density is minimal with respect to n p /n, which would lead to the constraint [59], where c 4 ≡ (2b 2 f 2 π m 2 K )/n 0 ∼ 49 MeV. After determining the ground state orientation, let us examine the consequences on the axion potential. The pseudo-NGB potential, after reintroducing the fluctuations we are mainly interested in, namely the neutral mesons π 0 , η and the axion, is given by withM a given in eq. (4.1) and  We consider densities in the region n < n c ≡ n 0 /b 1 ≈ 2.8n 0 (45 MeV/σ πN ) for the corresponding values of σ πN . The effect of kaon condensation is to eventually increase the neutral pion mass, while the axion becomes lighter and even massless at some density below n c . Above this density the axion field is therefore unstable around a = 0 and axion condensation is expected to occur.
in particular their decrease with density, ζū u ≈ ζd d ≈ 1 − b 1 (n/n 0 ). The three mass eigenstates, corresponding to mixtures of π 3 , η and a, have the following masses in the isospin symmetric limit ∆m = 0 and at leading order in θ andm/m s , 13 (m 2 π 0 ) n,θ ≈ m 2 π ζū u 1 + 1 8 (m 2 a ) n,θ ≈ (m 2 a ) 0 ζū u 1 − with m π and m η the masses in vacuum, respectively eq. (2.12) and m 2 η = −4m s qq 0 /3f 2 π neglecting O(m u,d /m s ) terms. Note that only the neutral pion mass depends on the charge chemical potential, and that such dependence enters along with kaon condensation. The effect of a non-vanishing θ on (m 2 π 0 ) n,θ therefore depends on the relative size of m s /m ≈ 27 and 2µ 2 /(m 2 π ζū u ), which enter with opposite signs. Note in this regard that while we use the leading order result for the quark condensate ratio ζū u , we did not perform an expansion in density. This is because, as we advanced at the beginning of this section and as explicitly shown in figure 2, when kaon condensation sets is we have n/n 0 > 1. Then, when m s /m > 2µ 2 /m 2 π ζū u , the coefficient of θ is negative and π 0 becomes lighter as kaon condensation sets in. Such a decrease could potentially lead to an instability and CP violation in the neutral sector. However, in the opposite case, when m s /m < 2µ 2 /(m 2 π ζū u ), which occurs at larger densities where ζū u is small and µ/m π large (see figure 2), the 13 The calculation of the axion mass is simplified by choosing a particular θ-dependent Qa matrix which removes the tree-level mixing between the axion and π0 and η, see appendix C for more details. coefficient of θ is positive and π 0 becomes heavier in the kaon-condensed phase. As opposed to the π 0 , the η mass depends only weakly on θ. In figure 3 we plot the numerical result for (m 2 π 0 ) n,θ as a function of density for fixed values of proton fraction, using the numerical results for θ(n) and µ(n) displayed in figure 2. We find that the µ 2 contribution leads to an increase in the mass of the neutral pion, similar to the effect of pion condensation in the simplified case discussed in section 3.2. Finally, the axion mass is independent of µ and decreases with the size of the kaon condensate. Interestingly, the negative coefficient of θ 2 is enhanced as density increases, since then ζū u becomes smaller. As shown in figure 4, this behavior eventually results in axion condensation at large densities, yet before the quark condensate vanishes. In this phase CP is thus spontaneously broken in the neutral sector.
Lastly, we note that one should be wary of the fact that for the quark condensates we included only density corrections at linear order and disregarded higher order corrections. As discussed at the end of the previous section, for densities n ∼ n c , these corrections are in fact important. In this regard, we would like to stress the fact that while our results might not be trustable at the quantitative level, this does not necessarily make an axion-condensed phase less likely. First, let us note that qualitatively we expect (m 2 a ) n,θ to decrease with n even when considering higher-order corrections to ζū u (n). Second, to further support our claim, let us consider the limit of maximal kaon condensation, i.e. θ → π/2, where If one ignores the density dependence of the condensates by taking b 1 = 0, this result is consistent with the naive expectation of a vanishing axion mass for cos θ → 0, since (m 2 a ) θ ∝ Tr[Σ 0 + Σ † 0 ] ∝ cos θ, as we showed in section 3.2. However, as discussed in section 4.1, a background of protons and neutrons makes ζū u < ζs s ≈ 1, in other words JHEP07(2020)221 b 1 > 0, this implies a negative axion mass for large kaon condensates, where it becomes energetically favorable for the axion field to develop a non-vanishing expectation value.

Axion couplings
Let us now turn our attention to the couplings of the axion to QCD matter at finite density. These couplings are of special importance, as they are a crucial input in the astrophysical axion bounds [14], in particular regarding supernovae and neutron star cooling, see e.g. [15,16,66,67] for recent works on the subject.
In vacuum. These coupling have been precisely calculated including RGE effects [25] and more recently at next-to-next-to-leading order in ChPT [68,69]. We review here how the couplings to protons and neutrons are derived in ChPT, following closely the discussion in [25]. The relevant part of the low-energy effective Lagrangian for the isospin doublet N = (p, n) T in the non-relativistic approximation is given by where we omitted mass terms proportional to the axion-dependent quark mass matrix M a , since there are no linear interactions coming from them if CP is preserved, see the discussion on CP violation at the end of this section. We also omitted higher-order terms in the expansion in (spatial) momenta, k/Λ χ ∼ k/M B . v µ is the velocity of the non-relativistic fields, which satisfy v µ γ µ N = N , while S µ is the spin operator, Here we have introduced the derivative ∇ µ , which contains the external (isospin) axial and vector currents, i.e. ∇ µ ξ = ∂ µ ξ−i(V µ −A µ )ξ and ∇ µ ξ † = ∂ µ ξ † − i(V µ + A µ )ξ † . 14 Finally,û i µ is associated to the (isospin) axial scalar current,û i µ = 2Â i µ , 15 where the index i = (u + d, s) runs over iso-scalar quark combinations.
Because of its UV couplings to the axial currents made out of quarks, eq. (2.9c), the axion enters eq. (4.40) as an external axial current, with components in both the iso-vector and iso-scalar directions. Therefore, with the inclusion of the axion, one finds with c IR ± ≡ (c IR u ±c IR d )/2 and where we have writtenû µ explicitly as a two-component vector. The couplings c IR u,d,s are related to the UV couplings by (4.43) 14 We have also introduced the Goldstone matrix ξ, because when the QCD orientation is trivial, ξL = ξ † R = ξ. Besides, note that in our convention uµ is twice that used in [25], following the standard in the ChPT literature, see e.g. [70] (however, in this reference the roles of ξ and ξ † are interchanged with respect to ours). 15 In case the η was light, thenû i µ ∝ ∂µη .

JHEP07(2020)221
where the matrix C qq accounts for renormalization group evolution (RGE) [25], see [25] for a detailed derivation. We recall from section 2 that the c 0 q are the UV model-dependent couplings of the axion to SM axial quark currents, while the matrix Q a was introduced in order to remove the aGG term. To further eliminate all axion-pion mixing, we chose a particular rotation matrix Q * a , which at zero density (denoted by the "0" subscript) is given by The coefficients g A and g i 0 in eq. (4.40) are given by linear combinations of hadronic matrix elements encoding the contribution of a quark q to the spin operator of the proton, The axial-vector coupling g A have been precisely measured in β-decay experiments, while g ud 0 and ∆s have been be extracted from lattice calculations, where in both cases isospin breaking effects are neglected. Their values at zero density are [25,71]   With this information, the couplings of the axion to nucleons, ∂ µ a f a (c pp S µ p + c nn S µ n) (4.49) can be finally extracted from the effective Lagrangian in eq. (4.40), This leads for instance to the accurate determination of the model-independent axion couplings, i.e. those of the KSVZ or hadronic axion [23,24], for which c 0 q = 0, We note in particular that the axion coupling to neutrons in vacuum is suppressed with respect to the naive O(1) expectation due to an accidental cancelation between z = m u /m d ≈ 1/2 and the ratio of matrix elements in vacuum ∆u/∆d = (g ud /(c n ) 0 ≈ 1.5. As we will be showing in the following, finite density corrections also spoil the cancellation, leading in fact to a much larger effect.
In-medium mixing angles. The mixing angles with the neutral pions change at finite baryon density due to the change in the values of the quark condensates, as discussed in section 4.1. The Q * a matrix that diagonalizes such mixings becomes therefore density dependent, where we defined Using eq. (4.12) for the quark condensates at linear order in density, we find where b 1,2,3 have been defined in terms of sigma terms in eq. (4.13). Note that the deviation of Z from unity is small, being proportional to the anomalously small coefficient b 2 , while the effects of W will be suppressed by m u,d /m s . However, similar to the RGE effects discussed above, the small density effect from Z gets amplified due to the cancellation structure of (c n ) KSVZ 0 , while no large enhancement is expected in c p . Indeed, one finds In-medium matrix elements. The hadronic matrix elements ∆u, ∆d and ∆s are also density dependent quantities in medium. In particular, the combination g A ≡ ∆u − ∆d, which fixes the coupling of the pions to nucleons, is known to get quenched inside nucleons [72]. This was observed from the reduced rates for β-decay in various large nuclei [73], suggesting that in medium (g A ) n 0 /2 /(g A ) 0 ≈ 0.75, with n 0 /2 being the typical baryon density around the surface of such large nuclei.

JHEP07(2020)221
The in-medium change in the effective axial coupling can be derived from the following higher-order terms in the non-relativistic baryon chiral Lagrangian [74,75] The density dependence of g A was originally calculated in [76] and used recently to explain the apparent discrepancy in β-decay rates in large nuclei [77]. 17 It is given, at leading order in n/Λ χ f 2 π (recall Λ χ ∼ M B ), by We identify two types of corrections. The terms proportional to I(m π /k f ) arise from (the resummation of) pion-exchange contributions originating in the operators in eq. (4.59), where k f = (3π 2 n/2) 1/3 ≈ (270 MeV) (n/n 0 ) 1/3 and we took the limit of vanishing momentum carried by the pion (there is little variation if instead the Fermi gas average value p 2 π = 6k 2 f /5 is taken). The contribution proportional to c D comes instead from the contact term in eq.
Similar to g A , higher-order operators in ChPT will give rise to a density dependent modification of g ud 0 ≡ ∆u + ∆d and ∆s,

66)
16 Note that in the literature these terms usually appear with dimensionful coefficients c3 and c4. Here we normalized them to the cutoff of ChPT, ci ≡ĉiΛχ. 17 This change in gA does not only affect the axion but also neutrino dynamics in supernovae (see e.g. [78] for a review on neutrinos in supernovae), which would be interesting to explore.

JHEP07(2020)221
from where one could derive, analogously to g A , the density corrections from pion exchange or contact terms. We parametrize our ignorance about the density dependence of the axialscalar coupling g ud 0 by defining κ, and, in light of eq. (4.64), we will consider the two benchmarks, κ = ±0.3, leading to in-medium quenching, as (g A ) n , or enhancement. Besides, while we could use a similar parametrization for ∆s, its contribution to the axion couplings to nucleons is already subleading, since (∆s) 0 = O(10 −2 ), thus we will neglect it in the following. We stress that it is certainly important to properly compute the density corrections to g ud 0 and ∆s, a challenging task given the expected uncertainties that would be associated with the determination of the coefficients in eqs. (4.65) and (4.66). We wish to point out however that we find no reason for the approximate relation ∆u/∆d ≈ −1/z ≈ −2 to hold even if all the relevant finite density corrections are taken into account.
Combining both the finite density effects discussed above, let us write the density dependent axion-nucleon couplings as the obvious generlization of eq. (4.51), where we recall that, since [(Q * a ) 0 ] q /[(Q * a ) n ] q ∼ b 2 (n/n 0 ) = O(10 −2 ) for q = u, d, the main effect of a baryonic background comes via (g A ) n (as well as (g ud 0 ) n ), a change that affects any type of axion (something fully encoded in the coefficients (c IR ± ) n ). To highlight the main point of this analysis, namely that the couplings of the axion to nucleons significantly change at finite density, we plot in figure 5 the ratio of the modelindependent but density-dependent axion couplings to nucleons (including RGE), normalized to the vacuum values, as a function of n/n 0 . As argued above, for such a hadronic axion the finite density effects are most striking, because the accidental suppression of the in-vacuum axion-neutron coupling is gone.
In the left panel, where we take κ positive, the O(1) modification of (g A ) n and (g ud 0 ) n translates into an O(10) enhancement of (c n ) KSVZ n at nuclear-saturation densities. In case of a negative κ, the ratio ∆u/∆d remains similar to its n = 0 value, and the accidental cancellation persists even after including in-medium effects leading only to O(1) modification of c n , although with large uncertainties. For the coupling to protons we find the opposite behaviour: for κ > 0 the increase in (g ud 0 ) n compensates for the decrease in (g A ) n , such that the coupling is almost unchanged at saturation density. On the other hand, for κ < 0, both (g 0 ) n and (g ud 0 ) n decrease, which leads to an O(1) decrease of c p .
In kaon-and axion-condensed phases. Let us briefly comment on the couplings of the axion when this acquires a non-trivial background value, which is the most interesting see e.g. [79,80]. While the effects of these couplings certainly deserve further investigation, in particular in the context of dense systems such as neutron stars, we only wish to point out that they are proportional to the quark masses, (4.72) since in the chiral limit the expectation value of the axion is not physical. 18

CFL phase
In this section we make a jump to asymptotically large baryon densities, or equivalently large quark chemical potentials, µ q Λ χ (µ q ≡ µ u = µ d = µ s ). At such high densities, two quark around the highly energetic Fermi surface interact weakly via tree-level gluon exchange, interactions which can be effectively parametrized below the Fermi momentum by 4-Fermi operators. Such operators, when in the attractive color3 channel, become relevant for back-to-back scattering as one approaches the Fermi surface (see e.g. [81,82]), leading to the formation of diquark pairs and ultimately to color superconductivity [31][32][33][34], see also [83] for a comprehensive review. Such a diquark pairing manifests itself in the form of a diquark condensate q L Cq L , which leads to the color-flavor-locked symmetry breaking pattern

JHEP07(2020)221
The condensates are given by [83] q ia L Cq jb L = abc ijk ∆3 L kc + ∆ 6 L ij,ab 3 We parametrize the low-energy fluctuations of the ground state, i.e. the NGBs associated with the symmetry breaking pattern (5.1) as The η and H are the NGBs associated with the spontaneous breaking of U(1) A and U(1) B , respectively. The NGBs associated with the breaking of color have been removed, since they are "eaten" by the gluons (unitary gauge). The rest of the NGBs, formally equivalent to those of the standard QCD chiral Lagrangian, are contained in the ξ L,R matrices, which are used to construct the following linearly-transforming color-neutral Goldstone matrix similarly as we did in vacuum, µ q = 0, see section 2. JHEP07(2020)221

Kinetic terms
The kinetic terms in the CFL phase are given by We recall that the introduction of chemical potential breaks Lorentz symmetry down to spatial rotations, and the low-energy excitations, even if massless, propagate sub-luminally. These velocities, as well as the decay constants, can be calculated by matching the UV microscopic theory [84] to the effective low-energy theory [85,86], 19 The Σ field gets a dynamically induced chemical potential due to the non-vanishing quark masses [88] Note that even if we choose a basis in which the axion enters the CFL effective Lagrangian via an axion-dependent quark mass matrix M a , as in eq. (2.9), it will not appear in such an effective chemical potential, since we can restrict ourselves to diagonal Q a matrices. In any case, for the analysis of the axion potential in the CFL phase, it will be more convenient to work in a basis where the axion is coupled to gluons, since a perturbative instanton expansion exists, being the gluons heavy and weakly coupled.

Mass terms
Given the spurionic transformation properties of the quark mass M in eq. (2.7), the leading order terms preserving the global symmetries in (5.1) are Note that this potential respects U(1) A and it is generated perturbatively. Contrary to QCD in vacuum, the axial symmetry thus dictates that the leading order terms in the scalar potential are O(M 2 ). Using eqs. (5.6)-(5.11) one finds [85] 20 (5.20) 19 As mentioned above, the gluons, with electric and magnetic masses m 2 D = g 2 s f 2 π and m 2 M = v 2 ϕ m 2 D respectively, are heavy and integrated out [87]. 20 The first term in eq. (5.19) can also be written as −2A1∆

JHEP07(2020)221
The coefficients can be computed by appropriately matching to the UV theory [85,86,89], Importantly, these two coefficients enter with opposite signs in eq. (5.19). This is due to the fact that while the color3 channel is attractive and lowers the total energy of the system, the color 6 channel is repulsive and increases it. As a result, one finds that ∆ 6 = 0 at the classical level. However, since ∆ 6 L,R does not break any additional symmetries compared to ∆3 L,R , there is nothing preventing it from being generated at the quantum level in the presence of a non-vanishing ∆ 3 . Indeed, a perturbative calculation yields [90] where α s = g 2 s /4π. ∆ 3 itself can be calculated using the so-called gap equation, in particular in the CFL phase with N f = 3 [83] The reason for considering the contribution of the condensate ∆ 6 to the potential, even though it is suppressed with respect to ∆ 3 , comes from the hierarchy in the quark masses. Indeed, one finds e.g. contributions from both condensates to the masses of the mesons, of order [35], The similarity of these two contributions, along with the fact that the coefficients of the respective operators (in eq. (5.21)) come with opposite signs, can lead to non-trivial vacuum structures, as we review below.

Non-perturbative terms
Instantons explicitly break the U(1) axial symmetry of QCD, also in the CFL phase [91]. At leading order in the gap parameters, one finds the following term generated via a single t'-Hooft vertex The coefficient A 3 can be calculated reliably at large chemical potentials due to the screening of gluons for instantons of size ρ 1/µ q 1/Λ QCD , where Λ QCD ≈ 250 MeV is the QCD scale parameter,

JHEP07(2020)221
with c = 0.155 [83,91]. Given that the operator in eq. (5.25) matches the leading term in the meson potential of the chiral Lagrangian at zero density, eq. (2.10), its coefficient can be mapped to the value of the standard quark condensate in the CFL phase where we set the chemical potential to a value expected to be realized in the core of a NS, noting that α s is to be evaluated at the scale µ q and that ∆ 3 depends on both α s and µ q . Due to the limited reliability of the perturbative result at such chemical potentials (see the discussion below), as well as to the strong dependence on Λ QCD , it is clear that one cannot make a robust prediction regarding the value of the quark condensate at realistic densities, yet a strong suppression of qq remains the most plausible outcome.
A higher-order operator that contributes to the mass of the η in the chiral limit appears at the two-instanton level, Note this term matches the would-be leading potential for the standard η in vacuum, see footnote 4. Before moving to the discussion of the axion potential in the CFL phase, let us note that the matching procedure by which the coefficients of the effective CFL Lagrangian are extracted from the microscopic theory relies on perturbative calculations that have been found to be under control for g s 0.8 [92]. Such a small coupling corresponds to very high quark chemical potentials, µ q 10 8 MeV, five orders of magnitude higher than those expected at the cores of dense NSs, where µ q ∼ 500 MeV. Still, a quantitative but more importantly a qualitative understanding of the CFL phase and of the corresponding axion potential provides a solid ground from which to extrapolate to lower chemical potentials and thus to realistic densities. In fact, the qualitative features and basic symmetry structure of the CFL phase should hold down to µ ∼ m 2 s /∆ 3 ≈ 180 MeV (50 MeV/∆ 3 ) [83].

Axion potential
In view of the previous discussion, in the following we examine the different axion potentials that arise by considering different hierarchies between the coefficients of the CFL operators.
Non-perturbative dominance. In this case we assume that the non-perturbative contributions to the potential dominate over the mass terms, that is V CFL 1-inst. , V CFL

2-inst.
V CFL M . Nevertheless, we still consider there exists a weak-coupling expansion, in the sense that the one-instanton contribution dominates over the two-instanton one, that is After a field redefinition η → η − (f η /4f a )a, this is found to be the same as in the vacuum chiral Lagrangian with a light η , which is minimized at the trivial vacuum, η = a = 0 in particular. The axion mass can be calculated by integrating out the mesons as we did in zero density. The details of this derivation can be found in appendix A. We find that the axion mass is suppressed with respect to its vacuum value by (5.31) Perturbative dominance. Let us now consider the hierarchy V CFL M V CFL

1-inst.
V CFL 2-inst. . One should first note that if the instanton terms are set to zero, the axion is massless, as can be immediately seen by using the basis defined by Q a = 0 in eq. (2.9); we use such a basis in the following. We write the potential in terms of the variables 4η /f η ≡ α , a/f a ≡ β , (5.32) and use the ansatz [35] where the angles ϕ and θ correspond to the expectation values η/ √ 3f π and K 0 /f π , respectively. In this basis the meson potential is given by  21 We should note at this point that we did not include the neutral pion π0 in our analysis because the corresponding first term in eq. (5.34a), which destabilizes the potential at the origin of field space for K 0 , vanishes for π0.

JHEP07(2020)221
The solution ϕ = α implies that the minimization of the NLO potential V CFL pert.,NLO with respect to α is found at Eqs. (5.36) and (5.37) match the results of [35]. Lastly, the instanton potential is minimized with respect to β at cos β = Sign [cos(3α)] . (5.38) Therefore, one finds that the axion is aligned with the η , such that while the axion mass, neglecting mixing with η and normalized to its vacuum value, is given by where we evaluated A 3 in eq. (5.26) at µ q = 1 GeV, Λ QCD = 250 MeV and α s = π. We therefore find that the axion can develop a non-vanishing expectation value in the CFL phase also, when the kaon condensate is large. Up to uncertainties associated with the value of A 3 , the axion is significantly lighter than in vacuum.

Axion sourcing observables
We briefly discuss in this section the potentially observable consequences of a non-vanishing axion condensate in NSs, where the largest baryonic densities among the stars are found. We defer to future work a more in-depth analysis of the corresponding phenomenology [30], as well as the study of the implications of the change in the axion-nucleon couplings with density, the latter particularly relevant for supernovae and NS cooling. For simplicity, let us consider the following toy model, namely a stepwise radiusdependent axion potential where m π and f π are the vacuum values and we have fixed the constants such that in the decoupling limit f a → ∞, the potential vanishes. The potential grossly captures the effect of matter on the axion potential, i.e. at a critical radius r c , which is of the order of the NS radius R, the axion field gets destabilized and the minimum of the potential is located at a/f a = π. The field equation can be solved numerically and one finds the intuitive result based on energy conservation, i.e. for the axion to get sourced the gain in potential energy needs to be enough to compensate for the gradient energy that comes with the change in In figure 6 we depict the typical field configurations of the axion sourcing, highlighting in grey the possible observable implication, to be discussed in turn below.

Free (vacuum) energy
The first potentially observable implication is associated with the shift in potential energy density inside the NS, as a result of the axion sourcing. This effect is independent of the field configuration outside the core of the NS, namely it is independent of m 2 a out . Such an energy density shift can be of considerable size compared to the energy density inside a NS, ρ 0 ≈ m n n 0 ≈ (190 MeV) 4 . Indeed, if the axion is sourced at relatively low baryon densities, as in the kaon condensed phase (section 4.2), one expects ∆V ∼ m 2 π f 2 π , which is indeed not significantly below ρ 0 or the energy change due to kaon condensation, of O(m 2 K f 2 π ). Instead, if axion sourcing happens in the CFL phase (section 5), this effect is expected to be suppressed by a few orders of magnitude, see eq. (5.40), and therefore likely negligible.

JHEP07(2020)221
A NS with a core of "vacuum energy" was considered as a generic scenario in [28,29], in the context of exotic QCD phases. It was found that the energy shift inside the NS leads to a significant change in the mass-radius relation of NSs, as well as to changes of the so-called chirp mass, one of main the gravitational wave observables of compact binary mergers. One could then, in a similar fashion, consider axion condensation as a possible source for such a potential energy shift.

Axion brane
Another interesting implication of axion condensation concerns the generation of a brane of energy density inside the NS. In particular, if the condition is met, the exponentially suppression of the axion field outside the NS is very rapid, a change that contributes to the energy density of the system in the form of a gradient energy One can think of such an abrupt change in field values as a localized (spherical) brane of energy density. The effect of such a brane has in fact not been considered in previous studies of NS structure. Such a brane would appear as an effective discontinuity in the temporal and radial components of the metric, which because of Einstein's field equations (known as the Tolman-Oppenheimer-Volkoff equations), imply a discontinuity in the pressure and enclosed mass of the star.

Axion-EM conversion
Next we consider the interplay between the EM fields of rotating NSs (i.e. pulsars), which are the strongest found in the Universe, and the axion, in particular when such that the sourced axion field is still non-negligible in the close surroundings of the NS. The axion and the classical EM fields form a coupled system, as seen from the generalized form of the Maxwell equations where the last line is the axion equation of motion. The interplay between the axion and the EM field of pulsars has been actively investigated before, see e.g. [93][94][95][96] for recent works on the subject, although the effects we consider here, associated with a large classical axion JHEP07(2020)221 field configuration also sourced by the NS, are novel. Assuming the conventional rotating dipole model, one finds that at the surface of the NS with Ω the angular velocity of the NS. Even with such large EM fields, we may still neglect the effects of the axion-photon coupling on the axion dynamics, since where we used (m 2 a ) out ∼ m 2 π (f π /f a ) 2 , a(R) ∼ f a and g aγγ ∼ α EM /f a . While we can safely assume that the back-reaction of the EM fields on the axion is negligible, it is also important to note that the value of a(R) decreases exponentially outside the NS, see eq. (6.3), such that one could well imagine a situation where the effect of the g aγγ (E · B) term is in fact comparable to the axion mass term. In this case the back-reaction of the EM fields would have to be taken into account, which is beyond the scope of this work.
We thus treat the axion field as a rigid source of additional EM fields, which can be simply estimated as thus leading to an O(1) enhancement around the surface of the NS. We note that since this correction is large, one could be concerned about whether the system can be treated perturbatively. This is in fact the case, since the higher order terms scale like This means that, apart from the leading O(α EM B * ) correction, further contributions are subleading. An additional sensitive observable is the dipole radiation output P that is responsible for the spin-down of rotating NSs. In this case we find that ∆P/P ∼ α 2 EM 1, namely there is no appreciable addition to the radiated energy due to the axion field. follows from the requirement of the axion being actually sourced, eq. (6.2). Therefore we expect the hierarchy (m a ) −1 in R (m a ) −1 out to arise only in non-standard scenarios, such as the one considered in [27]. If that is indeed the case, the long tails of the axion field configuration lead to a long range force between the NSs, generated by the Yukawa-like potential where Q eff = 4πf a R plays the role of the effective charge. This could lead to a deformation of the merger wave-form predicted by general relativity in case of NS with opposite-sign charges. A more dramatic effect would be found in the case of a repulsive force from same-sign charges, since in this case at some critical distance the axion force would dominate gravity, which could lead to halt in the merger process [27]. The presence of the axion field can also lead to an additional mechanism of energy loss in NS mergers, in the form of the scalar equivalent of Larmor radiation [97].

Conclusions and outlook
In this paper we have shown how the properties of the QCD axion change with baryon chemical potential. Based on the non-relativistic baryon chiral Lagrangian, reliable up to (slightly above) nuclear saturation densities, we have found that the axion gets lighter in medium, albeit the reduction is below an order of magnitude. In contrast, we have found that density corrections have a large impact on the axion couplings to nucleons, in particular the model-independent (a.k.a. KSVZ) coupling to neutrons is enhanced by up to one order of magnitude, thus becoming, in a baryonic background at nuclear saturation density, of the same order as the coupling to protons. We have also derived the conditions under which the axion acquires a non-trivial average value as a result of kaon condensation, an hypothetical yet motivated scenario for nuclear matter at densities above nuclear saturation.
Since reliable predictions for QCD matter at baryon densities as those found in the cores of neutron stars are currently unavailable, in order to delimitate the properties of the QCD axion not only from low but from high densities as well, we have analyzed the axion potential in the color-flavor-locked phase of QCD. We have again derived the conditions under which axion condensation takes place, and found that the axion mass is generically several orders of magnitude smaller than in vacuum.
The main goal of this paper was to lay the groundwork for future phenomenological investigations of axion physics in dense systems. This is the reason why we only brushed over some of the corresponding experimental implications, focussing in particular on the sourcing of the axion. In this context, we are currently studying in detail the implications of such an axion condensate contributing to the mass of neutron stars, as well as the idea of a brane separating the inner axion core from the rest of the neutron star [30], which are interesting idealizations of the richer dynamics of axions within these astrophysical objects. Another very compelling avenue concerns the interplay of the axion with the large electromagnetic fields generated by pulsars. Importantly, these are only a subset of the JHEP07(2020)221 range of exciting consequences our analysis has uncovered. In particular, the significance of the finite density corrections on the axion couplings calls for a proper assessment of how all such type of effects (as well as finite temperature) affect the predictions of the cooling rates of supernovae and proto-neutron stars, of special significance for the QCD axion.
Moreover, we believe our work will be of relevance for the study of other well-motivated new physics scenarios at finite density. This is certainly the case for those realizations of the so-called relaxion that rely on QCD dynamics -on the dependence of the quark masses on the electroweak scale -to generate a dynamical landscape for the Higgs fields. The investigation of this and similar density-dependent landscapes and its interplay with neutron stars is something that we leave for a future publication [40].

JHEP07(2020)221
Next we similarly integrate η out Finally, we integrate out η The potential is minimized around a = π 0 = η = η = 0 and we find the following axion mass Once we have diagonalized the mass matrix, one could be concerned with the effect of the O(f π /f a ) kinetic mixing which is generically induced by the derivative couplings of the axion. Let us explictly show that this does not affect the leading order results for the axion mass. Our starting point, without loss of generality, is the following Lagrangian and ξ ≡ f π /f a is our expansion parameter. Let start by performing the orthogonal rotation R 1 in the meson subspace, such that We rewrite the Lagrangian in this basis (A.9) We diagonalize and canonically normalize the lower 2 × 2 block in the second term by rotating and rescaling the fields

JHEP07(2020)221
T can be expanded T = 1 + 1 2 |d|ξσ 3 + O(ξ 2 ). At leading order T = 1 and the mass matrix can be re-diagonalized by performing the inverse orthogonal rotation R −1 2 , bringing it back to the diagonal form of eq. (A.7). One concludes that the axion mass does not receive any leading order correction due to the kinetic mixing.

(B.2)
It should be understood that the L and R operators act on different indices and therefore commute. We define the following rotated generators 3) The broken and unbroken generators are given by which transforms as Σ → LΣR † . Following standard notation, (B.8)

JHEP07(2020)221
We introduce the (θ-rotated) baryon octet as linearly-transforming fields,B θ L,R , where we use the θ-superscript because the finite-density backgrounds we consider consist of a non-vanishing ensemble of the standard (non-rotated) baryons, given bŷ with the usual parameterization The Lagrangian in this basis is given by where we recall the quark mass matrix spurion transforms as M → LM R † , we dropped some terms at the same order in derivatives (acting on the Σ matrices) that are irrelevant for our discussion, and we included leptons.

B.1 Adding chemical potential
We add chemical potentials for the three mutually commuting abelian symmetries associated with neutron and proton numbers and electric charge (from here on we neglect the other baryons). Under U(1) n,p , ψ → e iα ψ for ψ = n, p respectively, while under U(1) EM electromagnetism,  In this basis there are no non-derivative interactions of the mesons with the baryons from the mass terms in eq. (B.15). Besides, in complete analogy to eq. (B.10), the standard (non-rotated) baryons are given by In terms of such fields, which we recall make up the finite-density background, the baryon Lagrangian is given by Finally, we recall that at zero temperature all the states with E(p) = p 2 + m 2 ψ < µ ψ are occupied, such that with g ψ counting the internal degrees of freedom, e.g. g ψ = 2 for a fermion. In section 4.