Breakdown of chiral perturbation theory for the axion hot dark matter bound

We show that the commonly adopted hot dark matter (HDM) bound on the axion mass $m_a \lesssim$ 1 eV is not reliable, since it is obtained by extrapolating the chiral expansion in a region where the effective field theory breaks down. This is explicitly shown via the calculation of the axion-pion thermalization rate at the next-to-leading order in chiral perturbation theory. We finally advocate a strategy for a sound extraction of the axion HDM bound via lattice QCD techniques.

Introduction. The axion originally emerged as a low-energy remnant of the Peccei Quinn solution to the strong CP problem [1][2][3][4], but it also unavoidably contributes to the energy density of the Universe. There are two qualitatively different populations of relic axions, a non-thermal one comprising cold dark matter (DM) [5][6][7][8], and a thermal axion population [9] which, while still relativistic, would behave as extra dark radiation. Such hot dark matter (HDM) component contributes to the effective number of extra relativistic degrees of freedom [10] ∆N eff 4/7 (43/[4g S (T D )]) 4/3 , with g S (T D ) the number of entropy degrees of freedom at the axion decoupling temperature, T D . The value of ∆N eff is constrained by cosmic microwave background (CMB) experiments, such as the Planck satellite [11,12], while planned CMB Stage 4 (CMB-S4) experiments [13] will provide an observable window on the axion interactions.
There are several processes that can keep the axion in thermal equilibrium with the Standard Model (SM) thermal bath. From the standpoint of the axion solution to the strong CP problem, an unavoidable process arises from the model-independent coupling to gluons, αs 8π a fa GG. 1 For T D 1 GeV thermal axion production proceeds via its scatterings with gluons in the quark-gluon plasma [19,20], while for T D 1 GeV processes involving pions and nucleons must be considered [21][22][23]. The latter, have the advantage of occurring very late in the thermal history, so that it is unlikely that the corresponding population of thermal axions could be diluted by inflation. The transition between the two regimes depends on the strength of the axion interactions set by f a or, equivalently, by m a 5.7 × (10 6 GeV/f a ) eV, and it encompasses the range m a ∈ [0.01, 0.1] eV (with heavier axions leading to lower decoupling temperatures). Although the transition region cannot be precisely determined due to the complications of the quark-hadron phase transition, for heavier axions approaching the eV scale the main thermalization 1 Other thermalization channels arise from model-dependent axion couplings to photons [9], SM quarks [14][15][16][17] and leptons [18].
channel is aπ ↔ ππ [22,23], with T D 200 MeV. In this regime, scatterings off nucleons are subdominant because of the exponential suppression in their number density.
The highest attainable axion mass from cosmological constraints on extra relativistic degrees of freedom, also known as HDM bound, translate into m a 1 eV [24]. Based on a leading-order (LO) axion-pion chiral effective field theory (EFT) analysis of the axion-pion thermalization rate [22,23], the axion HDM bound has been reconsidered in Refs. [25][26][27][28][29][30][31][32][33], also in correlation with relic neutrinos. The most recent update [33] quotes a 95% CL bound that ranges from m a 0.2 eV to 1 eV, depending on the used data set and assumed cosmological model. Although the axion mass range relevant for the HDM bound is in generic tension with astrophysical constraints, the latter can be tamed in several respects. 2 It is the purpose of this Letter to revisit the axion HDM bound in the context of the next-to-LO (NLO) axion-pion chiral EFT. This is motivated by the simple observation that the mean energy of pions (axions) in a heat bath of T 100 MeV is E ≡ ρ/n 350 MeV (270 MeV), thus questioning the validity of the chiral expansion for the scattering process aπ ↔ ππ. The latter is expected to fail for √ s ∼ E π + E a 500 MeV, corresponding to temperatures well below that of QCD deconfinement, which was estimated to be T c = 154 ± 9 MeV in Ref. [43], see also [44,45].
In this work, we provide for the first time the formulation of the full axion-pion Lagrangian at NLO, including also derivative axion couplings to the pionic current (previous NLO studies only considered non-derivative axion-pion interactions [46,47]), and paying special attention to the issue of the axion-2 Tree-level axion couplings to electrons are absent in KSVZ models [34,35], thus relaxing the constraints from Red Giants and White Dwarfs. The axion coupling to photons, constrained by Horizontal Branch stars evolution, can be accidentally suppressed in certain KSVZ-like models [36][37][38]. Finally, the SN1987A bound on the axion coupling to nucleons can be considered less robust both from the astrophysical and experimental point of view [39][40][41][42].
pion mixing. Next, we perform a NLO calculation of the aπ ↔ ππ thermalization rate (that can be cast as an expansion in T /f π , with f π 92 MeV) and show that the NLO correction saturates half of the LO contribution for T χ 62 MeV. The latter can be considered as the maximal temperature above which the chiral description breaks down for the process under consideration. On the other hand, the region from T χ up to T c , where chiral perturbation theory cannot be applied, turns out to be crucial for the extraction of the HDM bound and for assessing the sensitivity of future CMB experiments.
We conclude with a proposal for extracting the axion-pion thermalization rate via a direct Lattice QCD calculation, in analogy to the well-studied case of π-π scattering.
Axion-pion scattering at LO. The construction of the LO axion-pion Lagrangian was discussed long ago in Refs. [48,49]. We recall here its basic ingredients (see also [22,50]), in view of the extension at NLO. Defining the pion Goldstone matrix U = e iπ A σ A /fπ , with f π 92 MeV, π A and σ A (A = 1, 2, 3) denoting respectively the real pion fields and the Pauli matrices, the LO axion-pion interactions stem from where χ a = 2B 0 M a , in terms of the quark condensate B 0 and the 'axion-dressed' quark mass matrix M a = e i a 2fa Qa M q e i a 2fa Qa , with M q = diag (m u , m d ) and Tr Q a = 1. The latter condition ensures that the axion field is transferred from the operator αs 8π a fa GG to the phase of the quark mass matrix, via the quark axial field redefinition q → exp(iγ 5 a 2fa Q a )q. In the following, we set Q a = M −1 q /Tr M −1 q , so that terms linear in a (including a-π 0 mass mixing) drop out from the first term in Eq. (1). Hence, in this basis, the only linear axion interaction is the derivative one with the conserved SU(2) A pion current. The latter reads at LO while the derivative axion coupling in Eq. (1) is where the first term arises from the axial quark rotation that removed the aGG operator and the second one originates from the model-dependent coefficient [34,35], while c 0 u = 1 3 cos 2 β and c 0 d = 1 3 sin 2 β in the DFSZ model [51,52], with tan β the ratio between the vacuum expectation values of two Higgs doublets. Expanding the pion matrix in Eq. (1) one obtains At the LO in the diagonalization of the a-π 0 term is obtained by shifting a → a − π 0 and π 0 → π 0 + O( 3 )a, where we used the fact that m a /m π = O( ). Hence, as long as we are interested in effects that are linear in a and neglect O( 3 ) corrections, the axionpion interactions in Eq. (3) are already in the basis with canonical propagators. For temperatures below the QCD phase transition, the main processes relevant for the axion thermalization rate are a(p 1 )π 0 (p 2 ) → π + (p 3 )π − (p 4 ), whose amplitude at LO reads with s = (p 1 + p 2 ) 2 , together with the crossed channels aπ − → π 0 π − and aπ + → π + π 0 . The amplitudes of the latter are obtained by replacing Taking equal masses for the neutral and charged pions, one finds the squared matrix element (summed over the three channels above) [23] Axion-pion scattering at NLO. To compute the axion thermalization process beyond LO we need to consider the one-loop amplitudes from the LO Lagrangian in Eq. (1) as well as the tree-level amplitudes stemming from the NLO axion-pion Lagrangian, both contributing to O(p 4 ) in the chiral expansion. The NLO interactions include the derivative coupling of the axion to the NLO axial current, which has been computed here for the first time.
We stick to the expression of the NLO chiral Lagrangian given in Ref. [53] (see for example Appendix D in [54] for trace notation), which, considering only two flavours, depends on 10 low-energy constants (LECs) 1 , 2 , . . . , 7 , h 1 , h 2 , h 3 . The axion field has been included in the phase of the quark mass matrix, as described after Eq. (1). Note that since we are interested in 2 → 2 scattering processes, we can neglect the O(p 4 ) Wess-Zumino-Witten term [55,56] since it contains operators with an odd number of bosons.
To compute the axial current J A µ at NLO, we promote the ordinary derivative to a covariant one, defined as D µ U = ∂ µ U − ir µ U + iU l µ , with r µ = r A µ σ A /2 and l µ = l A µ σ A /2 external fields which can be used to include electromagnetic or weak effects. The left and right SU(2) currents are obtained by differentiating the NLO Lagrangian with respect to l A µ and r A µ , respectively. Taking the R − L combination and switching off the external fields, the NLO axial current reads where curly brackets indicate anti-commutators. New a-π 0 mixings arise at NLO, both at tree level from the NLO Lagrangian and at one loop from L LO a-π . These mixings are explicitly taken into account in the Lehmann-Symanzik-Zimmermann (LSZ) reduction formula [57] (focussing e.g. on the aπ 0 → π + π − channel) where the index i runs over the external particles, Z a (Z π ) is the wave-function renormalization of the axion (pion) field and the full 4-point Green's function is given by The first term is the amputated 4-point function, multiplied by the 2-point functions of the external legs with the axion mass to zero. Working with LO diagonal propagators, the 2-point amplitude for the a-π 0 system reads P ij = diag (p 2 , p 2 − m 2 π ) − Σ ij , where Σ ij encodes NLO corrections including mixings. The 2-point Green's function G ij = (−iP) −1 ij is hence Plugging Eq. (9) and (10) into the LSZ formula for the scattering amplitude and neglecting O(1/f a ) 2 terms, one finds (with Z a = 1, Z π = 1 + Σ ππ (m 2 π ) and primes indicating derivatives with respect to p 2 ) where the G's are evaluated at the physical masses of the external particles. The one-loop amplitudes have been computed in dimensional regularization.
To carry out the renormalization procedure in the (modified) MS scheme, we define the scale independent parameters i as [53] i = with R = 2 d−4 − log(4π) + γ E − 1, in order to cancel the divergent terms (in the limit d = 4) with a suitable choice of the γ i . Eventually, only the terms proportional to 1,2,7 contribute to the NLO amplitude, which is renormalized for γ 1 = 1/3, γ 2 = 2/3 and γ 7 = 0. The latter coincide with the values obtained in Ref. [53] for the standard chiral theory without the axion.
The renormalized NLO amplitude for the aπ 0 → π + π − process (and its crossed channels) is given in Eq. (16) of the Supplementary Material. We have also checked that the same analytical result is obtained via a direct NLO diagonalization of the a and π 0 propagators, without employing the LSZ formalism with off-diagonal propagators. For consistency, we will only consider the interference between the LO and NLO terms in the squared matrix elements, , since the NLO squared correction is of the same order of the NNLO-LO interference, which we neglect.
Breakdown of the chiral expansion at finite temperature. The crucial quantity that is needed to extract the HDM bound is the axion decoupling temperature, T D , obtained via the freeze-out condition (following the same criterium as in [23]) Here, H(T ) = 4π 3 g (T )/45 T 2 /m pl is the Hubble rate (assuming a radiation dominated Universe) in terms of the Planck mass m pl = 1.22 × 10 19 GeV and the effective number of relativistic degrees of freedom, g (T ), while Γ a is the axion thermalization rate entering the Boltzmann equation where n eq a = (ζ 3 /π 2 )T 3 and f i = 1/(e Ei/T − 1). In the following, we will set the model-dependent axion couplings c 0 u, d = 0 (cf. Eq. (4)), to comply with the standard setup considered in the literature [22,23,[25][26][27][28][29][30][31][32][33] (see [58] for an exception). Moreover, we will neglect thermal corrections to the scattering matrix element, since those are small for T m π [59][60][61]. By integrating numerically the phase space in Eq. (14) we find (see Eq. (17)  step, or [62] for a slightly different approach)  [47], m u /m d = 0.50 (2) [64], f π = 92.1(8) MeV [24] and m π = 137 MeV (corresponding to the average neutral/charged pion mass). The h-functions are normalized to h LO (0) = h NLO (0) = 1 and they are plotted in Fig. 4 of the Supplementary Material. We have checked that h LO reproduces the result of Ref. [23] within percent accuracy. It should be noted that Eq. (15) is meaningful only for m π /T 1, since at higher temperatures above T c pions are deconfined and the axion thermalization rate should be computed from the interactions with a quark-gluon plasma. Nevertheless, we are interested in extrapolating the behaviour of Eq. (15) from the low-temperature regime, where the chiral approach is reliable.
In Fig. 1 we compare the LO and NLO rates contributing to Γ a = Γ LO a + Γ NLO a . In particular, the |Γ NLO a /Γ LO a | ratio does not depend on f a . Requiring as a loose criterium that the NLO correction is less than 50% of the LO one, yields T χ 62 MeV as the maximal temperature at which the chiral description of the thermalization rate can be reliably extended. Fig. 2 shows instead the extraction of the decoupling temperature (defined via Eq. (13)) for two reference values of the axion mass (setting the strength of the axion coupling via f a ), namely m a = 1 eV and 0.1 eV. Assuming a standard analysis employing the LO axion thermalization rate [23], the former benchmark (1 eV) corresponds to the most conservative HDM bound [33], while the latter (0.1 eV) saturates the most stringent one [33] and also represents the While in the former case the decoupling temperature is at the boundary of validity of the chiral expansion, set by T χ 62 MeV, in the latter is well above it. Hence, the region where the chiral expansion fails, T D T χ , corresponds to m a 1.2 eV. Since m a 1.2 eV yields a too large contribution to ∆N eff incompatible with Planck data (cf. Fig. 3), this value can be regarded as the axion HDM bound that can be reliably extracted within chiral perturbation theory. However, in the relevant mass range m a ∈ [0.1, 1] eV the decoupling temperature and consequently the axion HDM bound cannot be reliably extracted within the chiral approach.
Note, finally, that in the presence of modeldependent axion couplings c 0 u,d 1 (as in some axion models [65]), the same decoupling temperature as in the c 0 u,d = 0 case is obtained for larger f a , thus shifting down the mass window relevant for the axion HDM bound.
Towards a reliable axion HDM bound. The failure of the chiral approach in the calculation of the axion-pion thermalization rate can be traced back to the fact that in a thermal bath with temperatures of the order of T 100 MeV the mean energy of pions is E π 350 MeV, so that π-π scatterings happen at center of mass energies above the validity of the 2-flavour chiral EFT. The latter can be related to the scale of tree-level unitarity violation of π-π scattering resulting in √ s √ 8πf π 460 MeV [66,67]. A possible strategy to extend the theoretical predictions at higher energies is to compute the relevant aπ → ππ amplitudes using lattice QCD simulations. To this end one may employ the standard techniques used to compute weak non-leptonic matrix elements [68,69] and π-π scattering amplitudes as a function of the energy at finite volume [70][71][72][73]. Although this approach has limitations with respect to the maximum attainable center of mass energy, we believe that it can be used to compute the amplitudes up to values of √ s ∼ 600 − 900 MeV or higher [74]. We conclude by stressing the importance of obtaining a reliable determination of the axion-pion thermalization rate, not only in view of the extraction of a notable bound in axion physics, but also in order to set definite targets for future CMB probes of the axion-pion coupling, which could represent a 'discovery channel' for the axion.